NASA / CR— 1998-207946 



The Marshall Engineering Thermosphere 
(MET) Model 

Volume I: Technical Description 

R.E. Smith 

Physitron, Inc., Huntsville, Alabama 


Prepared for Marshall Space Flight Center 
under Contract NAS8-38333 

National Aeronautics and 
Space Administration 


Marshall Space Flight Center 




EXECUTIVE SUMMARY 


Volume I of this report presents a technical description of the Marshall Engineering 
Thermosphere (MET) model atmosphere and a summary of its historical development. Vari- 
ous programs developed to augment the original capability of the model are discussed in de- 
tail. The report also describes each of the individual subroutines developed to enhance the 
model. Computer codes for these subroutines are contained in four appendices. 

Volume II contains a copy of each of the reference documents. 

If you have questions or comments or need to request additional information, contact 
the Chief, Electromagnetics and Aerospace Environments Branch, Systems Analysis and Inte- 
gration Laboratory, MSFC, AL 35812. 


iii 




TABLE OF CONTENTS 


Section Title Page 


I INTRODUCTION 1 

n SCOPE 1 

m DESCRIPTION OF THE NEUTRAL THERMOSPHERE 1 

A. Variations 1 

B . Solar Activity Parameters 4 

IV HISTORICAL DEVELOPMENT 6 

A. Jacchia 1964 Model 6 

B. Jacchia 1970 Model 7 

C. Jacchia 1971 Model 7 

D. Marshall Engineering Thermosphere (MET) Model 7 

V BASIC MET MODEL 9 

A. Input Parameters for Real-Time or After-the-Fact Applications 10 

B. Input Parameters for Future Time Applications 10 

C. Output Parameters for All Applications 1 1 

VI BASIC MET MODEL WITH INTERNAL GRAVITY WAVES 1 1 

A. Optional Wave Models 12 

VII STATISTICAL ANALYSIS MODE, MET-S AM 13 

A. Input Parameters 14 

B. Output Parameters 14 

Vffl REFERENCES I 5 

APPENDICES 

A PROGRAM LISTING. MARSHALL ENGINEERING 

THERMOSPHERE (MET) MODEL A1 

B USER'S SOFTWARE IMPLEMENTATION GUIDE AND 

PROGRAM LISTINGS A SIMULATION OF SMALL-SCALE 
THERMOSPHERIC DENSITY VARIATIONS FOR ENGINEERING 
APPLICATIONS 

(For orbital inclinations <_40°) 

C USER'S SOFTWARE IMPLEMENTATION GUIDE AND 

PROGRAM LISTINGS AN ENGINEERING MODEL FOR 
A SIMULATION OF SMALL-SCALE THERMOSPHERIC 
DENSITY VARIATIONS FOR ORBITAL INCLINATIONS 
GREATER THAN 40° Cl 


V 



TABLE OF CONTENTS (continued) 


Title Pag e 

D USER'S SOFTWARE IMPLEMENTATION GUIDE AND 
PROGRAM LISTINGS THE MET MODEL STATISTICAL 
ANALYSIS MODE (MET-SAM) D1 


TABLE OF ILLUSTRATIONS 


Figure Title Page 

1 Variation in Geomagnetic Index, a p , in 1985 5 

2 Variation in F 10 7 Flux in 1985 6 


vi 



I. 


INTRODUCTION 


The Marshall Engineering Thermosphere (MET) Model of the Earth's atmos- 
phere at orbital altitudes was developed over a long period of time. This document 
summarizes the primary activities that occurred during the development and matu- 
ration of the model and the method and time various model options should be ap- 
plied in space vehicle development programs. 


II. SCOPE 

The first volume of this document includes a description of the Earth's neutral 
orbital atmosphere, summary of the historical development of the model, and de- 
scriptions of various options that can be exercised with the model. Detailed descrip- 
tions of these options, when and how they can be used, and program listings are in- 
cluded in the appendices. Reference documents are contained in Volume II. 


III. DESCRIPTION OF THE NEUTRAL THERMOSPHERE 

The region of the Earth's atmosphere between about 90 and 500 kilometers al- 
titude is known as the thermosphere. The region above about 500 kilometers is 
known as the exosphere. The temperature in the lower thermosphere increases rap- 
idly with increasing altitude from a minimum at 90 kilometers, but eventually be- 
comes altitude-independent at upper thermospheric altitudes. This asymptotic tem- 
perature, known as the exospheric temperature, does not vary with height at any 
given time due to the extremely short thermal conduction time. However, the exo- 
spheric temperature does vary with time because of solar activity and other factors 
discussed below. 

State of the neutral thermosphere is most conveniently described in terms of a mean, 
with spatial and temporal variations about that mean. The neutral thermosphere is 
important for two reasons. First, even at its low density, it produces torques and 
drag on orbiting spacecraft. Second, the density-height profile of the atmosphere 
above 100 kilometers altitude modulates the flux of trapped radiation encountered at 
orbital altitudes. 

A. VARIATIONS 

1. Solar Activity 

Short wavelength solar electromagnetic radiation (EUV and UV) 
changes substantially with level of solar activity. Thus, thermospheric density is 



strongly dependent on the level of solar activity. An average 11 -year solar cycle 
variation exists as well as a 27-day variation in density related to the average 27-day 
solar rotation period. Variations, however, tend to be slightly longer than 27 days 
early in the solar cycle when active regions occur more frequently at higher solar lati- 
tudes and slightly shorter than 27 days later in the solar cycle when the active re- 
gions occur more frequently closer to the Sun's equator. Coronal holes and active 
longitudes also affect this average 27-day variation. Changes in the thermospheric 
density related to changes in level of solar (and geomagnetic) activity, e.g., flares, 
eruptions, coronal mass ejections (CMEs), and coronal holes (CHs), can begin almost 
instantaneously(minutes to hours), although more often a lag of a day or more oc- 
curs. 


2. Geomagnetic Activity 

Interaction of solar wind with the Earth's magnetosphere (referred to 
as geomagnetic activity) leads to a high latitude heat and momentum source for the 
thermospheric gases. Some of this heat and momentum is convected to low lati- 
tudes. Geomagnetic activity varies, usually has one peak in activity just prior to and 
another just after the peak activity of the solar cycle as defined by the 10.7-cm solar 
radio noise flux. Also, larger solar cycle peaks are associated with more intense 
geomagnetic activity. A seasonal variation of geomagnetic activity occurs with 
maxima in March (+ 1 month) and September (+ 1 month) each year. This variation 
is possibly related to the tilt of the Sun's rotational axis toward the Earth. 

3. Diurnal 

Rotation of the Earth induces a diurnal (24-hour period) variation 
(diurnal tide) in thermospheric temperature and density. Due to a lag in response of 
the thermosphere to the EUV heat source, density maximizes around 2 p.m. local 
solar time at orbital altitudes at a latitude approximately equal to the subsolar point. 
The lag decreases with decreasing altitude. Similarly, minimum density occurs be- 
tween 3 and 4 a.m. local solar time at about the same latitude in the opposite hemi- 
sphere. In lowest regions of the thermosphere (120 km and below), where character- 
istic thermal conduction time is on the order of a day or more, the diurnal variation is 
not a predominant effect. 

Harmonics of the diurnal tide are also induced in the Earth's atmosphere. A semi- 
diurnal tide (period of 12 hours) and a ter-diurnal tide (period of 8 hours) are impor- 
tant in the lower thermosphere (below about 160 km for the semi-diurnal tide and 
much lower for the ter-diurnal tide). Because of large damping effects of molecular 
viscosity, these diurnal harmonic tides are not important at orbital altitudes. 
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4. Semiannual 


This variation is believed to be a conduction mode of oscillation 
driven by a semiannual variation in Joule heating in the high latitude thermosphere 
(as a consequence of a semiannual variation in geomagnetic activity). The variation 
is latitudinally independent and is modified by compositional effects. Amplitude of 
the variation is height dependent and variable from year to year with a primary 
minimum in July, primary maximum in October, secondary minimum in January, 
followed by a secondary maximum in April. Magnitude and altitude dependence of 
the semiannual oscillation vary considerably from one solar cycle to another. This 
variation is important at orbital altitudes. 

5. Seasonal-latitudinal in Lower Thermosphere 

These seasonal-latitudinal variations are driven in the thermosphere 
by the dynamics of the lower atmosphere (mesosphere and below). Amplitude of the 
variation maximizes in the lower thermosphere between about 105 and 120 kilome- 
ters and diminishes to zero around 200 kilometers. Although the temperature oscil- 
lation amplitude is quite large, corresponding density oscillation amplitude is small. 
This variation is not important at orbital altitudes. 

6. Seasonal-latitudinal of Helium 

Satellite mass spectrometers measured a strong increase in helium 
above the winter pole. Over a year the helium number density varies by a factor of 
42 at 275 km, 12 at 400 km, and 3 or 4 above 500 km. Formation of the winter helium 
bulge is primarily due to effects of global scale winds that blow from the summer to 
the winter hemisphere. Amplitude of the bulge decreases with increasing levels of 
solar activity due to increased effectiveness of exospheric transport above 500 km 
carrying helium back to the summer hemisphere. Also, a very weak dependence ex- 
ists of helium bulge amplitude on magnitude of the lower thermospheric eddy dif- 
fusivity. 


7. Waves 

Atmospheric waves have been detected in temperature and density 
measurements throughout the atmosphere from the ground to at least 510 km. These 
fluctuations, caused by gravity waves, are so named because they are primarily oscil- 
lations of the neutral gas for which the restoring force is gravity. A thermospheric 
gravity wave produces a corresponding wave in the ionosphere known as a traveling 
ionospheric disturbance. 

Thermospheric gravity waves oscillate with periods typically of 30 minutes to several 
hours and have horizontal wavelengths from hundreds of kilometers to about 4000 
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km. Density amplitudes of the larger scale waves are greater at higher latitudes and 
smaller toward the equator. At approximately 200 km altitude typical values for 
density wave amplitudes are 15% at auroral latitudes and 5% at equatorial latitudes. 
Gravity wave amplitudes generally decrease at greater altitudes in the thermosphere 
due to dissipation by molecular processes. Larger scale waves survive to greater alti- 
tudes than the smaller. 

8. Winds 

At low latitudes (less than 28.5°) wind speeds range from 100 to 200 
meters per second. At high latitudes(greater than about 65 ) speeds can be 1500 me- 
ters per second or more. Rapid (minutes) changes in wind direction (to 180°), 
probably driven by gravity waves, have been observed by a satellite. 

B. SOLAR ACTIVITY PARAMETERS 

Various surrogate indices are used to quantify levels of solar activity. One is 
the 10.7-cm solar radio noise flux, designated F 107 . Although EUV radiation heats the 
atmosphere, this radiation cannot be measured on the ground. The F 107 can be meas- 
ured on the ground and correlates well with the EUV radiation. 

An index used as a measure of geomagnetic activity is the planetary geomagnetic ac- 
tivity index, a p , (or k p - essentially the logarithm of a p ). This index is based on mag- 
netic fluctuation data reported every 3 hours at 12 stations between geomagnetic lati- 
tudes 48 and 63 selected to provide good longitude coverage. Although high lati- 
tude ionospheric current fluctuations drive the magnetic field fluctuations observed 
at these stations, the magnetic field fluctuations do not drive the thermosphere. 

Thus, good correlations are not always found between observed density changes and 
the a p index. 

Figures 1 and 2 show the variability, during a period of low solar activity, of these 
two indices (a p and F 107 ). 
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Figure 1. Variation in Geomagnetic Index, a p/ in 1985. 
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Figure 2. Variation in F 107 Flux in 1985. 


IV. HISTORICAL DEVELOPMENT 

A. TACCHIA 1964 MODEL 

In the middle 1960s personnel at the NASA Marshall Space Flight Center 
(MSFC) responsible for predicting the orbital decay histories of satellites began 
studies to determine which atmospheric model when combined with the appropriate 
orbit propagation program would most accurately predict the observed decay histo- 
ries of 39 satellites. Since the observed decay histories were already available, a 
number of available models of the thermosphere were selected for use in appropriate 
computer programs. Because the satellites had decayed prior to the start of the 
study, actual values of proxy input parameters required by the models were used. 
These proxy parameters were representative of actual solar conditions that occurred 
during the decay periods of the satellites. The Jacchia 1964 model had the best per- 
formance statistically and, therefore, was selected for use by MSFC. 
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B. TACCHIA 1970 MODEL 


Almost coincidental with the completion of MSFC's study the Smithsonian As- 
trophysical Observatory (SAO) published the Jacchia 1970 model atmosphere (ref. 1). 

Due to the great similarity between the 1964 and 1970 models and confirmation from 
observed data of the existence of major differences, the decision was made to use the 
newer model. In the 1964 model the temperature induced density bulge remained on 
the equator all year while the 1970 model bulge followed the latitudinal excursions of 
the Sun. 

C. TACCHIA 1971 MODEL 

In 1971 SAO published the Jacchia 1971 model (ref. 2). Although this model 
had several new features confirmed by observational data, overall, the model was 
not as representative of the atmosphere as the 1970 model. 

D. MARSHALL ENGINEERING THERMOSPHERE (MET) MODEL 

An independent study by environmental personnel at MSFC showed that the 
1970 model could be improved by adding seasonal-latitudinal variations in the lower 
thermosphere density below 170 kilometers and adding seasonal-latitudinal varia- 
tions in helium from the 1971 model above 500 kilometers altitude while retaining 
the orbital decay prediction accuracy of the 1964 model. Subsequent minor fairing 
modifications eliminated step function increases in density spuriously introduced by 
these modifications. This was the first version of the MET model used at MSFC. 
Additional minor modifications made in the late 1980s corrected programming er- 
rors in earlier computer programs and made the complete program more under- 
standable and user friendly. For publications on programming changes and com- 
puter codes for the changed model see references 3, 4, and 5. 

1. Statistical Analysis Mode, MET-SAM 

This was the status quo for the MET model at the beginning of the 
Space Station Program when a meeting was held at the NASA Johnson Space Center 
(JSC) to establish design criteria for the development, testing, and operational phases 
of the Station. The complexities of the vehicle, requirement for a guaranteed 30-year 
orbital lifetime, and reboost strategy required to fulfill this requirement, plus the re- 
quirements of the various users, made this a most formidable task faced with very 
tight monetary restrictions. A complicated reboost strategy had to be developed to 
protect the astronauts from radiation at orbital altitudes and to ensure a 30-year life- 
time without outside assistance. A limit existed on the altitude attainable by the 
Space Shuttle as well as the number of flights the Shuttle could make for on-orbit 
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construction of the Station and for resupply of fuel and reboost energy. This strategy 
affected the size of the propellant tanks on-board Space Station and Shuttle to trans- 
port additional fuel to the Station. Adequate precautions also had to be taken to en- 
sure that the Space Station would not re-enter if the Shuttle missed a scheduled re- 
boost event. 

A critical factor in these deliberations at JSC was not only the representativeness of 
the atmospheric model but also the values of the input parameters which control 
model output. The atmospheric model required the 162-day mean values of 10.7-cm 
solar radio noise flux, previous day value of the 10.7-cm solar radio noise flux, and 3- 
hourly value of the geomagnetic index, K p or a p/ 6.7 hours prior to time of analysis 
application. For any future time application the only available values of these three 
parameters were predictions of the 13-month smoothed values of the 10.7-cm solar 
radio noise flux and the geomagnetic index, a p , based on statistics of the data during 
all previous solar cycles. Although statistics of separate distributions of each of the 
three parameters existed, nothing was available on probabilities of combined occur- 
rences of the three parameters. As a result of this void, the decision at the JSC meet- 
ing was to use different values of these two parameters, the 13-month smoothed val- 
ues of 10.7-cm solar radio noise flux and the a p during development, testing, and op- 
eration phases of the Space Station program. 

Early in the Station developmental phase, the need for a detailed investigation of the 
combined occurrence of these three parameters became evident. A study was begun 
to determine exactly how these three proxy input parameters should be combined 
for use in future time period activities. The result of this study is the Statistical 
Analysis Mode, MET-SAM, described in reference 6. 

2. Small Scale Thermosphere Simulations 

Studies of the Guidance, Navigation & Control (GN&C) Systems for 
Space Station revealed some of the more recently discovered short time period phe- 
nomena in the thermosphere (not included in the MET model) were critical to the 
GN&C systems performance. The existence of waves in the thermosphere was dis- 
covered in the OGO 6 satellite density measurements. When the MET model was 
used in GN&C studies involving a space vehicle having widely separated centers of 
gravity and pressure, the effects of these waves were very significant. This signifi- 
cance was further enhanced by operational constraints placed on the Space Station by 
users — notably the micro-gravity researchers' requirements for 30 consecutive days 
of gravity levels below 10‘ 6 gs. Since short period fluctuations could be the result of 
fluctuations in solar activity, the probability of occurrence of 30-day quiescent peri- 
ods was a concern. Environmental specialists decided to add to MET a wave model. 
An in-depth analysis was made of observed density measurements to construct this 
numerical model of the wave spectra expected to occur during all levels of solar ac- 
tivity. Since design activities progressed rapidly and answers were needed for on- 
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going engineering analyses before the observational data were completely analyzed, 
representative sets were selected from the 10.7-cm solar activity and a p databases for 
five different levels of solar activity that might exist during orbital operations of the 
Space Station. Data sets of both parameters were selected and tailored to reflect the 
95 percentile profile that might exist during selected 90-day periods (length of time 
scheduled between Space Shuttle reboost and refueling missions). Results of this 
analysis of the interim procedure to be used in the engineering analyses until the in- 
depth wave study is completed are contained in reference 7. 


3. Waves 

Available observational data from the AE-E satellite in a low inclination 
(19.7°) orbit were obtained from the National Space Sciences Data Center, GSFC, 
while programs were being written to accomplish the required analysis. The wave 
model described in reference 8 is a callable subroutine addition to the MET model. 

4. Waves2 

Almost at completion of this wave study and subsequent inclusion in 
the MET model, the announcement was made that the Space Station would be a joint 
effort between the USA and Russia and in a 51.6° inclination orbit. All previous 
analyses indicated most waves present in the thermosphere originated near the auro- 
ral ovals, propagated equatorward with phase fronts aligned with constant latitudes 
and with amplitudes that decreased as they propagated. Analyzed observational 
data from satellites in higher inclination orbits indicated a new wave model should 
be developed for inclusion in the MET model and applied to the International Space 
Station. The analysis and resulting callable subroutine for the MET model are de- 
scribed in reference 9. 


V. BASIC MET MODEL 

The basic MET model, ref. 3, is an empirical static diffusion model with coef- 
ficients obtained from satellite drag analyses. With the proper input parameters, 
specified below, the exospheric temperature can be calculated. In the original devel- 
opment phase of the model the prime objective was to model the total neutral mass 
density of the thermosphere by adjusting temperature profiles until agreement be- 
tween modeled and measured total densities was achieved. Agreement between 
modeled temperatures and temperatures measured on later missions was not always 
achieved. Thomson-scatter radar temperature measurements generally show that 
the diurnal temperature maximum lags the density maximum by a couple of hours, 
whereas in the MET model the temperature and the density maxima and minimal are 
in phase. 
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With exospheric temperature specified, the temperature can be calculated for any al- 
titude between 90 and 2500 km from an empirically determined temperature profile. 
The density for all points on the globe at 90-km altitude is assumed constant and 
mixing of atmospheric constituents prevails to 105 km. Between these two altitudes 
the mean molecular mass varies as a result of dissociation of molecular oxygen to 
atomic. At 120-km altitude the ratio of atomic oxygen to molecular is assumed to be 
1.5. Density between 90 and 105 km is calculated by integration of the barometric 
equation. For altitudes above 105 kilometers the diffusion equation for each of the 
individual species (O , O, N , He, and Ar) is integrated upward from the 105-km 
level. For hydrogen (H) integration of the diffusion equation proceeds upward from 
500-km altitude. Total mass density is calculated by summing the individual specie 
mass densities. 

The total mass density is modified further by the effects of seasonal-latitudinal den- 
sity variation of the lower thermosphere below 170-km altitude and seasonal- 
latitudinal variations of helium (He) above 500km. These two effects were incorpo- 
rated into the MET model using the equations developed by Jacchia for his 1971 
thermospheric model (ref. 2). 

The basic MET model includes all variations listed in Section III, with exception of 
the wind and wave variations. A single computer program exists for all applications 
of the basic MET model; however, two applications, listed below, require different 
values of the input parameters, also listed below. Appendix A is a listing of the 
FORTRAN source code for this basic MET model. A more detailed description of the 
model is in references 3, 4, and 5. 

A. INPUT PARAMETERS FOR REAL-TTME OR AFTF.R-THE-FACT 

APPLICATIONS 

1. Time: Year, month, day, hour, minute 

2. Position: Altitude, geographic latitude and longitude 

3. Solar Activity Parameters: Previous day's value of the 10.7-cm solar 
radio noise flux, centered value of the 10.7-cm solar radio noise flux averaged over 6 
solar rotations (162 days), a p or K p geomagnetic index 6 to 7 hours prior to time of 
application. 

B. INPUT PARAMETERS FOR FUTURE TIME APPLICATIONS 

1. Time: Year, month, day, hour, minute 

2. Position: Altitude, geographic latitude and longitude 

3. Solar Activity Parameters At present the only available data are 
predicted 13-month smoothed values of both the 10.7-cm solar radio noise flux and 
geomagnetic index, a p . 
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C. OUTPUT PARAMETERS FOR ALL APPLICATIONS 


1. Exospheric temperature 

2. Temperature 

3. N number density 

4. O number density 

5. O number density 

6. Ar number density 

7. He number density 

8. H number density 

9. Average molecular weight 

10. Total mass density 

11. Log mass density 

12. Total pressure 

13. Local gravity acceleration 

14. Ratio specific heats 

15. Pressure scale height 

16. Specific heat constant p 

17. Specific heat constant v 

The total mass density, temperature, and individual species number densities all 
have the same phase variation in the MET model. 


VI. BASIC MET MODEL WITH INTERNAL GRAVITY WAVES 

The shortest time period variations in the basic MET model occur every three 
hours as a result of the ap variations, giving a Nyquist period of six hours. Tides 
with periods of twelve and eight hours are known to exist at thermospheric altitudes 
and have small amplitudes. The smallest time period for all of these variations far 
exceeds the average orbital period for near-Earth orbiting vehicles. The only known 
variations at thermospheric altitudes with periods less than an orbital period are 
those associated with internal gravity waves. 

Atmospheric gravity waves exist from the troposphere to the top of the thermo- 
sphere. Their spatial and temporal scales are commensurate with their sources and 
both scales increase with increasing altitude. These waves can transport energy and 
momentum over large distances. Thermospheric gravity waves generally have spa- 
tial and temporal scales greater than 100 km and 20 minutes, respectively; however, 
waves with periods of several hours and wavelengths of several thousand kilometers 
have been inferred from observations. 
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Those waves, generated in situ at high latitudes in the thermosphere in response to 
geomagnetic activity, have larger scales and phase speeds. The phase speeds, as they 
propagate predominantly toward the equator, depend on the level of geomagnetic 
activity. 

A. OPTIONAL WAVE MODELS 

Optional additions to the basic MET model realistically model the observed 
wave spectra. Results of analyses of orbital altitude density measurements show to- 
tal mass density in the thermosphere is best represented by summing the basic MET 
model total mass density and a quasi-random wave component having an absolute 
amplitude proportional to the mean basic MET model density. 

1 . FOR VEHICLES WITH ORBITAL INCLINATIONS LESS THAN OR 
EQUAL TO 40° The FORTRAN source code is listed in Appendix B. See reference 8 
for additional details. 

a. Input Parameters 

1. All input parameters for the basic MET model 

2. 13-month smoothed values of the 10.7-cm solar radio noise flux 
and geomagnetic index, A p 

b. Output Parameters 

1. Amplitude of the wave expressed as a percentage of and to be 
added to the basic MET model total mass density 

2. FOR VEHICLES WITH ORBITAL INCLINATIONS GREATER THAN 
40 Source code is listed in Appendix C. See reference 9 for additional details. 

a. Input Parameters 

1. All input parameters for the basic MET model 

2. A daily A p value 

b. Output parameters 

1. Amplitude of the wave expressed as a percentage of and to be 
added to the basic MET model total mass density 
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VII. STATISTICAL ANALYSIS MODE, MET-SAM 


The Statistical Analysis Mode (SAM) of the MET model was developed for use 
primarily during the design, development, and testing phases of a space vehicle 
when only statistical estimates of the 13-month smoothed mean values of the three 
solar parameters are required as inputs. While the exospheric temperature in MET- 
SAM is obtained in a different manner than in the basic MET model, densities for the 
same exospheric temperature are identical in both. 

The MET-SAM provides answers to frequently asked questions during design and 
development of a space vehicle— what is the percentage of time that the recom- 
mended density value will be exceeded during the operational phase of the vehicle? 
Or, how confident are you that the recommended density value will not be exceeded 
more than a specified percentage of time? The answers are crucial in the design of 
the guidance and control capability, selection of the altitude at which the vehicles 
will orbit, amount of shielding necessary to protect the crew and vehicle from solar 
radiation and orbital debris, and reboost strategy for payloads put in low-Earth orbit 
by the Shuttle. MET-SAM is based on the premise that most applications during the 
development phase require detailed knowledge about maximum and minimum 
densities that will be encountered with limits on the magnitude of variations that oc- 
cur during monthly time periods. With a slight modification MET-SAM has the ca- 
pability to include the diurnal (daily) variation provided the input solar parameters 
remain constant during the orbit. 

Initially a new database was constructed for 1947 through 1991. Contents were the 
eight 3-hourly values of a p , daily value of 10.7-cm solar radio noise flux, and 162-day 
mean value of 10.7-cm solar radio noise flux centered on day of application refer- 
enced to the 13-month smoothed mean value of 10.7-cm solar radio noise flux cen- 
tered on month of application. The global minimum, mean, and maximum exo- 
spheric temperatures were calculated for every three-hour period from 1947 through 
1991 using algorithms in the MET model and appropriate solar activity input pa- 
rameters from this new database. These temperatures were sorted into five levels of 
solar activity, as defined by the 13-month smoothed value of the 10.7-cm solar radio 
noise flux. Cumulative percentage frequency (CPF) distributions were calculated for 
these temperatures and levels of solar activity. Results of these calculations are given 
in Tables 1 through 18, Appendix D. Range of temperature in each distribution is the 
result of various combinations of three solar activity parameters that are required in- 
put parameters for the basic MET model. All uncertainties in how the three parame- 
ters should be combined are included in the distributions in the tables based on how 
they occurred in the historical records. 

This procedure allows the design engineer to select the risk level he is willing to ac- 
cept by selecting the appropriate temperature and density from the CPF distribu- 
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tions. A complete description is in reference 6 and computer program listing in Ap- 
pendix D. 

To use the MET-SAM: 

A. INPUT PARAMETERS 

1. Altitude, in kilometers 

2. Appropriate temperature from CPF distributions tables 

B. OUTPUT PARAMETERS 

1. Exospheric temperature 

2. Temperature at altitude 

3. N number density 

4. O number density 

5. O number density 

6. Ar number density 

7. He number density 

8. H number density 

9. Average molecular weight 

10. Total mass density 

11. Log mass density 

12. Total pressure 

13. Local gravity acceleration 

14. Ratio specific heats 

15. Pressure scale height 

16. Specific heat constant p 

17. Specific heat constant v 
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£************ ****** ****** ******** ********** ********************* ************* 


c* * 

C* The Marshall Space Flight Center 

C* Marshall Engineering Thermosphere Model 

C* * 

c* * 


C* written by * 

C* * 

C* Mike Hickey * 

C* Universities Space Research Association * 

C* NASA / MSFC , ED44 * 

C* Tel. (205) 544-5692 * 

C* * 

C* This program is a driving program for the following subroutines 
C* * 

C* ATMOSPHERES * 

C* SOLSET * 

C* TIME * 

C* J70 * 

C* * 


C* The atmospheric model is a modified Jacchia 1970 model and is given in * 
C* the subroutine J70. All of the other subroutines were designed to * 

C* allow flexible use of this model so that various input parameters could * 

C* be varied within a driving program with very little software development.* 
C* Thus, for example, driving routines can be written quite easily to * 

C* facilitate the plotting of output as line or contour plots. Control is * 

C* achieved by setting the values of four switches in the driving program, * 

C* as described in subroutine ATMOSPHERES. * 

C* * 


** 


REAL*4 INDATA (12) , OUTDATA (12) , AUXDATA (5) 
CHARACTER* 1 SWITCH (4) 

C CALL LIB$INIT_TIMER 

C Set all switches to 'Y' so that only one particular calculation is performed 


SWITCH (1) = 'Y' 

SWITCH (2) = 'Y' 

SWITCH (3) = 'Y' 

SWITCH (4) = 'Y' 

CALL ATMOSPHERES ( INDATA, OUTDATA, AUXDATA, SWITCH ) 
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C Now type output data 


= OUTDATA (4),' /m3' 
= ', OUTDATA (5),' /m3' 
= OUTDATA (6),' /m3' 
= \ OUTDATA (7),’ /m3' 
= OUTDATA (8),' /m3' 


WRITE! 1,*)' All output in MKS units' 

WRITEU,*)" 

WRITE(1,*)' Exospheric temperature = OUTDATA (1),’ K' 
WRITE! 1,*)' Temperature = OUTDATA (2),' K' 

WRITE( 1 ,*)' N2 number density = OUTDATA (3),' /m3’ 
WRITE(1,*)' 02 number density 
WRITE(1,*)' O number density 
WRITE(1,*)' A number density 
WRITE( 1,*)' He number density 
WRITE! 1,*)' H number density 
WRITE! 1,*)' Average molecular wt. = OUTDATA (9) 

WRITE! 1,*)' Total mass density = OUTDATA (10),' kg/m3’ 
WRITE! 1,*)' Log 10 mass density = ', OUTDATA (1 1) 
WRITE! 1,*)' Total pressure = ', OUTDATA (12),' Pa' 

WRITE! 1,*)' Local grav. acceln. = AUXDATA (1),' m.sec-2' 
WRITE! 1 ,*)' Ratio specific heats = ', AUXDATA (2) 
WRITE(1,*)' Pressure scale-height = ', AUXDATA (3),' m' 
WRITE! 1,*)' Specific heat cons, p = ', AUXDATA (4),'m2.sec-2.K-l’ 
WRITE! 1,*)' Specific heat cons, v = ', AUXDATA (5),'m2.sec-2.K-r 
WRITE! 1,*)” 


CALL LIB$SHOW_TIMER 


STOP 

END 


SUBROUTINE ATMOSPHERES ( INDATA, OUTDATA, AUXDATA, SWITCH ) 


£**************************************************************************** 

** 

C* DESCRIPTION:- * 

Q* * 

c* * 

C* Calculate atmospheric data in single precision using subroutine J70 * 

C* and J70SUP. * 

c * * 

C* SUBROUTINES:- * 

c * * 

C* * 

C* TIME, SOLSET, GMC, J70 and J70SUP * 

* 

C* INPUT:- * 
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C* 

C* all single precision, either through 

C* subroutines or from main driver prog. 

C* * 

C* INDATA (1)- altitude 


(2) — latitude 

(3) — longitude 

(4) - year (yy) 

(5) — month (mm) 

(6) — day (dd) 

(7) — hour (hh) 

(8) — mins (mm) 


= Z 

= XLAT 
= XLNG 
= IYR 

= MN 

= IDA 
= IHR 

= MIN 


C* 

C* 

c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 

C* OUTDATA (1) -- exospheric temperature (K) 


(9) — geomagnetic index = IGEO_IND 

(10) -- solar radio noise flux= F10 

(11) — 1 62-day average F 1 0 = F 1 OB 

(12) — geomagnetic activity index = GI=AP 


OUTPUT- 


NOTE : All output in MKS units 


all single precision 


c* 

.. (2) — temperature at altitude Z 


* 

c* 

.. (3) — N2 number density (per meter-cubed) 


c* 

.. (4) — 02 number density ( 

) 

* 

c* 

.. (5) — O number density ( 

) 

* 

c* 

.. (6) — A number density ( 

) 

* 

c* 

.. (7) — He number density ( 

) 


c* 

.. (8) — H number density ( 

) 

* 

c* 

.. (9) — average molecular weight 

* 


c* 

.. (10)— total density 

* 


c* 

.. (11)— log 10 ( total density ) 

* 


c* 

c* 

.. (12)— total pressure ( Pa ) 

* 



C* AUXDATA ( 1) — gravitational acceleration ( m/s-s ) 


c* . 

. (2) — ratio of specific heats 

* 

C* . 

. (3) — pressure scale-height ( m ) 

* 

C* . 

. (4) — specific heat at constant pressure 

* 

c* . 

. (5) — specific heat at constant volume 

* 

c* 

* 


c* 

* 


c* 

COMMENTS:- 


c* 
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C* * 

C* SWITCH(l) — if Y(es), date and time are input from terminal through * 

C* subroutine TIME once only * 

C* SWITCH(2) -- if Y(es), solar/magnetic activity are input from terminal * 

C* through subroutine SOLSET once only 

C* SWITCH(3) -- if Y(es), only ONE altitude value is input from terminal * 
C* through main calling program * 

C* SWITCH(4) — if Y(es), only ONE latitude AND longitude are input from * 
C* terminal through main calling program * 

C* * 

C* ATMOSPHERES written by Mike Hickey ( USRA, NASA/ED44 ) 

C* Tel: (205) 544-5692 * 

C* January-April 1987 * 


^**************************************************************************** 


EXTERNAL TIME 
DIMENSION AUXDATA (5) 

INTEGER HR 

REALM LAT, LON, INDATA (12), OUTDATA (12) 
CHARACTER*1 SWITCH (4) 

PARAMETER PI = 3.14159265 


C 

C This next section is only executed on the first call to ATMOSPHERES 

DO WHILE ( CALL. EQ. 0.0 ) 

C SECTION A:- 
C 


IF ( SWITCH(l). EQ. 'Y' ) THEN 

CALL TIME ( IYR, MON, IDA, HR, MIN, SWITCH(l) ) 
INDATA (4) = FLO AT J (IYR) 

INDATA (5) = FLO AT J (MON) 

INDATA (6) = FLO AT J (IDA) 

INDATA (7) = FLOATJ (HR) 

INDATA (8) = FLOATJ (MIN) 

END IF 

C SECTION B:- 
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C 


IF ( SWITCH(2). EQ. 'Y' ) THEN 

CALL SOLSET ( IGEOJND, F10, F10B, GI, SWITCH(2) ) 
INDATA (9) = FLOATJ (IGEOJND) 

INDATA (10) = F10 
INDATA (11) = F10B 
INDATA (12) = GI 

END IF 

SECTION C:- 


IF ( SWITCH(3). EQ. 'Y’ ) THEN 

TYPE Input altitude, km' 
ACCEPT *, INDATA (1) 

Z = INDATA (1) 

END IF 

SECTION D:- 


IF ( SWITCH(4). EQ. 'Y' ) THEN 

TYPE Input latitude and longitude, degrees' 

ACCEPT *, ( INDATA(I), 1= 2,3 ) 

LAT = INDATA (2) 

LON = INDATA (3) 

RLT = INDATA (2) * PI / 180. ! geographic latitude, radians 
END IF 

CALL = 1.0 

END DO 

End of first executable section 


C The following depend on the values of the switches 
£*** + 

C* SECTION 1:- 

IF ( SWITCH(l). NE. 'Y' ) THEN 
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IYR = JINT ( INDATA (4) ) 

MON = JINT ( INDATA (5) ) 

IDA = JINT ( INDATA (6) ) 

HR = JINT ( INDATA (7) ) 

MIN = JINT ( INDATA (8) ) 

CALL TIME ( IYR, MON, IDA, HR, MIN, SWITCH(l) ) 

END IF 

£***** 

C* SECTION 2:- 

IF ( SWITCH(2). NE. Y' ) THEN 

IGEO _IND = JINT ( INDATA (9) ) 

F10 = INDATA (10) 

F10B = INDATA (11) 

GI = INDATA (12) 

CALL SOLSET ( IGEOJND, F10, F10B, GI, SWITCH(2) ) 
END IF 

C* SECTION 3:- 

IF ( SWITCH(3). NE. T' ) THEN 
Z = INDATA (1) 

END IF 

£***** 

C* SECTION 4:- 

IF ( SWITCH(4). NE. Y’ ) THEN 

LAT = INDATA (2) 

LON = INDATA (3) 

RLT = INDATA (2) * PI / 180. ! geographic latitude, radians 
END IF 

C All setting-up complete. 
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CALL J70 ( INDATA, OUTDATA ) 

CALL J70SUP ( Z, OUTDATA, AUXDATA ) 


RETURN 

ENTRY ATMOS_ENT ( DUMMY ) 

CALL = DUMMY 

RETURN 

END 


SUBROUTINE TIME ( IYR, MON, IDA, HR, MIN, SWITCH ) 


£******************** K*************************************************^****^ 


** 

c* * 

C* This subroutine sets up time of year and day * 

C* * 

C* INPUTS/OUTPUTS: * 

C* * 


C* IYR = year ( 2 digits ) * 

C* MON = month * 

C* IDA = day of month * 

C* HR = hour of day * 

C* MIN = minutes * 

C* * 

C* Written by Mike Hickey, USRA * 


c**************************************************************************** 

** 


DIMENSION IDAY ( 12) 
INTEGER HR 
CHARACTER*! SWITCH 


PARAMETER PI = 3.14159265 
DATA IDAY/ 31, 28, 31,30, 31,30, 31,31,30, 31,30, 31 / 

If SWITCH = Y(es) then input data and time from terminal 
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C. 


IF ( SWITCH.EQ.'Y'. OR. SWITCH.EQ.'y' ) THEN 

TYPE *, ' Input date and time of date? ( yy,mm,dd,hh,mm ) ' 
ACCEPT *, IYR, MON, IDA, HR, MIN 

END IF 


C 


IF ( JMOD (IYR, 4) .EQ. 0 ) THEN 

IF ( JMOD (IYR, 100) .NE. 0 ) ID AY ( 2 ) = 29 

ELSE 

IDAY ( 2 ) = 28 

END IF 
DAYTOT = 0.0 

DO 1 1= 1, 12 

DAYTOT = DAYTOT + FLOATJ ( IDAY ( I ) ) 

1 CONTINUE 

IF ( MON. GT. 1 ) THEN 


KE = MON - 1 
ID = 0 

DO 2 1= 1, KE 
ID = ID + IDAY (I) 

2 CONTINUE 


ID = ID + IDA 

DD = IDA 

END IF 


ELSE 


RETURN 

END 


SUBROUTINE SOLSET ( IGEOJND, F10, FI OB, GI, SWITCH ) 
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^^c*************************************************************************** 

*♦ 

c* * 

C* This subroutine simply calls for a setup of the solar-activty and auroral * 

C* activity indices. 

C* 

C* INPUTS/OUTPUTS: 

C* 

C* IGEO_IND = geomagnetic index 
C* F10 = solar radio noise flux 

C* FI OB = 162-day average FI 0 
C* GI = geomagnetic activity index 
C* 

C* Written by Mike Hickey, USRA 
** 



CHARACTER* 1 SWITCH 
IGEOJND = 2 


C 

C If SWITCH = Y(es) then input geomagnetic indices from terminal 

C 


IF ( SWITCH.EQ.'Y'. OR. SWITCH.EQ.’y’ ) THEN 

C TYPE *, ’ Input geomagnetic index ( 1-KP, 2-AP ) ' 

C ACCEPT *, IGEOJND 

TYPE *, ' Input solar radio noise flux ( F10 = 0-400 ) ' 
ACCEPT *, F10 

TYPE *, ' Input 162-day average F10 ( F10B = 0-250 ) ' 
ACCEPT*, FI 0B 

C IF ( IGEOJND . EQ. 2 ) THEN 

C 

C TYPE *, 1 Input geomagnetic activity index ( GI = 0-400 ) ' 
C 

C ELSE 

C 

C TYPE *, ' Input geomagnetic activity index ( GI = 0-9 ) ' 

C 

C END IF 

TYPE *,' Input AP index ( AP = 0 - 400 ) ' 
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ACCEPT *, GI 
END IF 


RETURN 

END 


SUBROUTINE J70SUP ( Z, OUTDATA, AUXDATA ) 


£**** ******* sit******)*!********* ******** ****** ***** ************** ****** ********* 
** 

c* * 

C* DESCRIPTION:- * 

C* * 

c* * 

C* J70SUP calculates auxiliary variables which are output in array * 

C* AUXDATA, given data input from J70 which are contained in array OUTDATA. * 

C* * 


c* 

INPUT DATA:- 

* 

c* 


* 

c* 


* 

c* 

Z — altitude (km) 

* 

c* 

T 77 . — temperature at altitude z = 

OUTDATA (2) 

c* 

— N2 number density = 

•• (3) 

c* 

- 02 .. .. = .. (4) 

* 

c* 

- 0 .. = .. (5) 

* 

c* 

- A .. = .. (6) 

* 

c* 

- He .. .. = .. (7) 

* 

c* 

- H .. = .. (8) 

* 

c* 

EM — average molecular weight 

= •• (9) 

C* DENS — total density = 

.. (10) * 

c* 

P — total pressure = 

(12) 

c* 


* 

c* 

OUTPUT DATA:- 

* 

c* 


* 

c* 


* 

c* 

G -- gravitational acceleration 

= AUXDATA (1) 

C* GAM — ratio of specific heats 

= AUXDATA (2) 

c* 

H — pressure scale-height 

= AUXDATA (3) 

c* 

CP -- specific heat at constant pressure = AUXDATA (4) 


All 
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C* CV -- specific heat at constant volume = AUXDATA (5) * 

C* * 

C* Written by Mike Hickey, USRA * 

C** ****** ****** ******** ****************************************************** 

** 


REAL*4 OUTDATA (12), AUXDATA (5), H 

G = 9.80665 / ( ( 1 . + Z / 6.356766E3 )**2 ) 

H = OUTDATA ( 1 2) / ( G * OUTDATA (10)) 

SUM1 = OUTDATA (3) + OUTDATA (4) 

SUM2 = 0.0 

DO 1 1 = 5,8 

SUM2 = SUM2 + OUTDATA (I) 

1 CONTINUE 

GAM = ( 1.4 * SUM1 + 1.67 * SUM2 ) / ( SUM1 + SUM2 ) 

CV = G * H / ( (GAM- 1.0) * OUTDATA (2) ) 

CP = GAM * CV 

AUXDATA (1) = G 
AUXDATA (2) = GAM 
AUXDATA (3) = H 
AUXDATA (4) = CP 
AUXDATA (5) = CV 

RETURN 

END 


SUBROUTINE J70 ( INDATA, OUTDATA ) 


£****** 

£** 

c** 

c** 

c** 

c** 

c** 

£** 

c** 


********************************************************************** 

** 

J70 developed from J70MM by ** 

Mike P. Hickey ** 

Universities Space Research Association ** 

at ** 

NASA / Marshall Space Flight Center, ED44, ** 

Huntsville, Alabama, 35812, USA. ** 

Tel. (205) 544-5692 ** 
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through the subroutine calling list 

** 

through the subroutine calling list 

** 




INPUT DATA: 


** 


Z - altitude 
XLAT — latitude 
XLNG -- longitude 
IYR - year (yy) 

MN — month (mm) 

IDA — day (dd) 

IHR - hour (hh) 

MIN — mins (mm) 

II — geomagnetic index 
F10 — solar radio noise flux = INDATA (10) 
F10B — 162-day average F10 = INDATA (1 1) 

GI — geomagnetic activity index = INDATA (12) 


= INDATA (1) 

= INDATA (2) 

= INDATA (3) 
= INDATA (4) 

= INDATA (5) 
= INDATA (6) 

= INDATA (7) 

= INDATA (8) 
= INDATA (9) 


** 


** 


** 


** 


** 

** 


C** 

C** INPUTS: 

C** 

C** OUTPUTS: 
c** 

c** 
c** 
c** 
c** 
c** 
c** 

Q ** 

C** 

c** 
c** 
c ** 
c** 
c** 
c** 

c** 
c** 

£** 
c** 
c** 

c** 

£** 
c** 
q** 

c** 

£** 

£** 

£** 
c** 

c** 

c** 
c** 

Q** 

£** 

Q** 

Q**************************************************************************** 
** 


OUTPUT DATA: 


** 


T— exospheric temperature = OUTDATA (1) 
TZZ— temperature at altitude Z = OUTDATA (2) 


A(l)-- N2 number density 
A(2)-- 02 number density 
A(3)~ O number density 
A(4)- A number density 
A(5)— He number density 
A(6)— H number density 
EM— average molecular weight = OUTDATA (9) 
DENS— total density = OUTDATA (10) 
DL— loglO ( total density ) = OUTDATA (1 1) 

P— total pressure = OUTDATA (12) 


= OUTDATA (3) 
= OUTDATA (4) 
= OUTDATA (5) 
= OUTDATA (6) 
= OUTDATA (7) 
= OUTDATA (8) 


** 

** 


** 


NB. Input through array 'INDATA' 
Output through array ’OUTDATA’ 




DIMENSION A ( 6 ) 

REAL*4 INDATA ( 12 ), OUTDATA ( 12 ) 

PARAMETER RGAS = 8.31432E3 ! J/kmol-K 

PARAMETER BFH = 440.0 


C Calcultions performed for only one latitude , one longitude 
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C and one altitude 
C 

C** Set parameters to INDATA values 
C 

Z = INDATA (1) 

XLAT = INDATA (2) 

XLNG = INDATA (3) 

IYR = JINT ( INDATA (4) ) +1900 
MN = JINT ( INDATA (5) ) 

IDA = JINT ( INDATA (6) ) 

IHR = JINT ( INDATA (7) ) 

MIN = JINT ( INDATA (8) ) 

II = JINT ( INDATA (9) ) 

F10 = INDATA (10) 

F10B = INDATA (11) 

GI = INDATA (12) 


CALL TME ( MN , IDA , IYR , IHR , MIN , XLAT , XLNG , SDA , 
SHA , DD , DY ) 

CALL TINF ( F10 , F10B , GI , XLAT , SDA , SHA , DY , II , TE ) 

CALL JAC ( Z , TE , TZ , A(l) , A(2) , A(3) , A(4) , A(5) , A(6) , 

EM , DENS , DL ) 


DENLG = 0. 

DUMMY = DL 
DEN = DL 

IF ( Z .LE. 170.) THEN 

CALL SLV ( DUMMY , Z , XLAT , DD ) 
DENLG = DUMMY 

END IF 


** 'Fair' helium number density between base fairing height ( BFH ) and 500 km 


IF ( Z. GE. 500. ) THEN 

CALL SLVH ( DEN , A(5) , XLAT , SDA ) 
DL = DEN 

ELSE IF ( Z ,GT. BFH ) THEN 
DHEL1 = A ( 5 ) 

DHEL2 = A ( 5 ) 

DLG1 = DL 
DLG2 = DL 
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CALL SLVH ( DLG2 , DHEL2 , XLAT , SDA ) 

IH = Z 

CALL FAIR5 ( DHEL1 , DHEL2 , DLG1 , DLG2 , IH , FDHEL , FDLG ) 
DL = FDLG 
A ( 5 ) = FDHEL 

END IF 


DL = DL + DENLG 

DENS = 10.**DL 

XLAT = XLAT * 57.29577951 


C Fill OUTDATA array 

OUTDATA (1) = TE 
OUTDATA (2) = TZ 

DO 801= 1,6 

OUTDATA (1+2) = 1.E6 * ( 10. ** A(I) ) 
80 CONTINUE 

OUTDATA (9) = EM 
OUTDATA (10) = DENS * 1000. 
OUTDATA (11) = DL 
P = OUTDATA (10) * RGAS * TZ / EM 
OUTDATA (12) = P 

RETURN 

END 


SUBROUTINE TME ( MN , IDA , IYR , IHR , MIN , XLAT , XLNG , 

SDA , SHA , DD , DY ) 

£**************************************************************************** 

** 

C** Subroutine 'TME' performs the calculations of the solar declination ** 

C** angle and solar hour angle. ** 

£** ** 

C** INPUTS: MN = month ** 

C** IDA = day ** 

C** IYR = year ** 

C** IHR = hour ** 

C** MIN = minute ** 

C** XMJD= mean Julian date ** 

C** XLAT= latitude ( input-geocentric latitude ) ** 
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C** XLNG= longitude ( input-geocentric longitude, -180,+ 180 ) 


C** ** 

C** OUTPUTS: SDA = solar declination angle (rad) ** 

C** SHA = solar hour angle (rad) ** 

C** DD = day number from 1 JAN. ** 

C** DY = DD / tropical year ** 

C** Modified by Mike Hickey, USRA ** 


£********************************************, It******************************* 
** 


DIMENSION IDAY(12) 


PARAMETER YEAR = 365.2422 
PARAMETER A1 = 99.6909833 
PARAMETER A2 = 36000.76892 
PARAMETER A3 = 0.00038708 
PARAMETER A4 = 0.250684477 
PARAMETER B 1 = 0.0172028 
PARAMETER B2 = 0.0335 
PARAMETER B3 = 1 .407 
PARAMETER PI = 3.14159265 
PARAMETER TPI = 6.2831853 1 
PARAMETER PI2 = 1.57079633 
PARAMETER PI32 = 4.71238898 
PARAMETER RAD_DEG = 0.017453293 
DATA IDAY / 3 1 ,28,3 1 ,30,3 1 ,30 ,31,31,30,31 ,30,3 1 / 

XLAT = XLAT / 57.2957795 1 
YR = IYR 

IF ( JMOD(IYR,4) .EQ. 0 ) THEN 

IF ( JMOD(IYR, 100) .NE. 0 ) IDAY(2) = 29 ! Century not a leap year 
ELSE 

IDAY(2) = 28 
END IF 

ID = 0 

IF ( MN. GT. 1 ) THEN 
DO 20 1= 1 , MN-1 

ID = ID + IDAY(I) 

20 CONTINUE 
END IF 

ID = ID + IDA 
DD = ID 
DY = DD/YEAR 
C 

C** Compute mean Julian date 
C 

XMJD = 2415020. + 365. * ( YR - 1900. ) + DD 
+ FLOATJ ( (IYR - 1901 )/4) 
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c 

C** Compute Greenwich mean time in minutes GMT 
C 

XHR = IHR 

XMIN = MIN 

GMT = 60 * XHR + XMIN 

FMJD = XMJD - 2435839. + GMT / 1440. 

C 

C** Compute Greenwich mean position - GP ( in rad ) 

C 

XJ = ( XMJD - 2415020.5 ) / ( 36525.0 ) 

GP = AMOD ( A1 + A2 * XJ + A3 * XJ * XJ + A4 * GMT , 360. ) 


C 

C** Compute right ascension point - RAP ( in rad ) 

C 

C** 1st convert geocentric longitude to deg longitude - west neg , + east 
C 

IF ( XLNG .GT. 1 80. ) XLNG = XLNG - 360. 

RAP = AMOD ( GP + XLNG , 360. ) 


C 

C** Compute celestial longitude - XLS ( in rad ) - - zero to 2PI 
C 

Y1 = B1 * FMJD 

Y2 = 0.0 17202* (FMJD -3.) 

XLS = AMOD ( Y1 + B2 * SIN(Y2) - B3 , TPI ) 

C 

C** Compute solar declination angle - SDA ( in rad ) 

C 

B4 = RAD_DEG * ( 23.4523 - 0.013 * XJ ) 

SDA = ASIN ( SIN ( XLS ) * SIN ( B4 ) ) 

C 

C** Compute right ascension of Sun - RAS ( in rad ) - - zero to 2PI 

C 

c** These next few lines do not appear in NASA CR- 179359 or NASA CR -?????? 
C** They are added here to ensure that that argument of ASIN stays bounded 
C** between -1 and +1 , which could otherwise be effected by roundoff error. 

ARG = TAN ( SDA ) / TAN ( B4 ) 

IF( ARG .GT. 1.0) ARG = 1.0 
IF (ARG .LT. -1.) ARG = -1.0 
RAS = ASIN ( ARG ) 

C 

C** Put RAS in same quadrant as XLS 
C 

RAS = ABS ( RAS ) 

TEMP = ABS ( XLS ) 


A17 



MARSHALL ENGINEERING THERMOSPHERE MODEL PROGRAM LISTING 

(PAGE 17 OF 28) 


IF ( TEMP.LE.PI .AND. TEMP.GT.PI2 ) THEN 
RAS = PI - RAS 

ELSE IF ( TEMP.LE.PI32 .AND. TEMP.GT.PI ) THEN 
RAS = PI + RAS 

ELSE IF ( TEMP.GT.PI32 ) THEN 
RAS = TPI - RAS 

END IF 

IF ( XLS. LT. 0. ) RAS = -RAS 


C 

C** Compute solar hour angle - SHA ( in deg ) - - 
C 

SHA = RAP * RAD_DEG - RAS 
IF ( SHA.GT.PI ) SHA = SHA - TPI 
IF ( SHA.LT.-PI ) SHA = SHA + TPI 


RETURN 

END 


SUBROUTINE TINF ( F10 , FI OB , GI, XLAT, SDA , SHA , DY , II ,TE ) 


c** 

c** 

c** 

c** 

c** 


** 

** 


** 




C** Subroutine 'TINF calculates the exospheric temperature according to 
C** L. Jacchia S AO 3 1 3, 1 970 ** 

£** ** 

F10 = solar radio noise flux ( x E-22 Watts / m2 ) ** 

F10B= 162-day average F10 
GI = geomagnetic activity index 
LAT 35 geographic latitude at perigee ( in rad ) 

SDA = solar declination angle ( in rad ) 

C** SHA = solar hour angle ** 

C** dy = D / Y ( day number / tropical year ) ; 1 
C** II = geomagnetic equation index ( 1— GI=KP , 2— GI=AP ) 

C** RE = diurnal factor KP t F10B, AVG 

** 

CONSTANTS — C = solar activity variation 
— BETA , etc = diurnal variation 
— D = geomagnetic variation 
— E = semiannual variation 

** 


c ** 

C** 

c ** 

c** 

c** 

c** 


** 

** 

** 


C** Modified by Mike Hickey, USRA 
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^^^^^^Jlc^ + Jt:*^********************************************************** ******* 

** 


PARAMETER PI = 3.14159265 

PARAMETER TPI = 6.283 18531 
PARAMETER XM = 2.5 

PARAMETER XNN = 3.0 
C 

C** Ci are solar activity variation variables 
C 

PARAMETER Cl = 383.0 

PARAMETER C2 = 3.32 
PARAMETER C3= 1.80 
C 

C** Di are geomagnetic variation variables 
C 

PARAMETER Dl = 28.0 
PARAMETER D2 = 0.03 
PARAMETER D3 = 1.0 
PARAMETER D4 = 100.0 
PARAMETER D5 = -0.08 
C 

C** Ei are semiannual variation variables 
C 

PARAMETER El =2.41 

PARAMETER E2 = 0.349 
PARAMETER E3 = 0.206 
PARAMETER E4 = 6.2831853 
PARAMETER E5 = 3.9531708 
PARAMETER E6 = 12.5663706 
PARAMETER E7 = 4.3214352 
PARAMETER E8 = 0. 1 145 

PARAMETER E9 = 0.5 
PARAMETER E10 = 6.2831853 
PARAMETER El 1 = 5.9742620 
PARAMETER E12 = 2.16 

PARAMETER BETA = -0.6457718 
PARAMETER GAMMA = 0.7504916 
PARAMETER P = 0. 1047 198 
PARAMETER RE = 0.31 
C 

C** solar activity variation 
C 

TC = Cl + C2 * F10B + C3 * ( F10 - F10B ) 
C 

C** diurnal variation 
C 

ETA = 0.5 * ABS ( XLAT - SDA ) 
THETA = 0.5 * ABS ( XLAT + SDA ) 
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TAU = SHA + BETA + P * SIN ( SHA + GAMMA ) 

IF ( TAU. GT. PI ) TAU = TAU - TPI 
IF ( TAU. LT.-PI ) TAU = TAU + TPI 

A1 = ( SIN ( THETA ) )**XM 
A2 = ( COS ( ETA ) )**XM 
A3 = ( COS ( TAU / 2. ) )**XNN 
B 1 = 1 .0 + RE * A1 
B2 = ( A2-A1 )/B! 

TV = B 1 * ( 1 . + RE * B2 * A3 ) 

TL = TC * TV 
C 

C** geomagnetic variation 
C 

IF ( Il.EQ.l ) THEN 

TG = D1 * GI + D2 * EXP(GI) 

ELSE 

TG = D3 * GI + D4 * ( 1 - EXP ( D5 * GI ) ) 

END IF 
C 

C** semiannual variation 
C 

G3 = 0.5 * ( 1.0 + SIN ( E10 * DY + El 1 ) ) 

G3 = G3 ** E12 

TAU1 = DY + E8 * ( G3 - E9 ) 

Gl = E2 + E3 * ( SIN ( E4 * TAU1 + E5 ) ) 

G2 = SIN ( E6 * TAU 1 + E7 ) 

TS = El + F10B * G 1 * G2 
C 

C** exospheric temperature 
C 

TE = TL + TG + TS 

RETURN 

END 


SUBROUTINE JAC ( Z , T , TZ , AN , A02 , AO , AA , AHE , AH , EM , 

DENS , DL ) 

^Hc****^**************************^*******^************************** ********* 
** 

( 2 ** ** 

C** Subroutine 'JAC' calculates the temperature TZ , the total density DENS ** 
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C** and its logarithm DL, the mean molecular weight EM, the individual ** 

C** specie number densities for N, 02, O, A, HE and H ( each preceded with ** 

€** an 'A' ) at altitude Z given the exospheric temperature T. ** 

C** This subroutine uses the subroutine 'GAUSS' and the function ** 

C** subprograms 'TEMP' and 'MOL_WT'. ** 

* * 

C** Rewritten by Mike Hickey, USRA ** 

£**************************************************************************** 

** 


DIMENSION ALPHA(6) , EI(6) , DI(6) , DIT(6) 

REAL*4 MOL_WT 

PARAMETER AV = 6.02257E23 
PARAMETER QN =.78110 
PARAMETER Q02 = .20955 
PARAMETER QA =.009343 
PARAMETER QHE = 1.289E-05 
PARAMETER RGAS = 8.31432 
PARAMETER PI =3.141 59265 
PARAMETER TO =183. 

C GRAVITY ( ALTITUDE ) = 9.80665/((l. + ALTITUDE/6.356766E3)**2) 


DATA ALPHA / 0.0 , 0.0 , 0.0 , 0.0 , -.380 , 0.0 / 

DATA El / 28.0134 , 31.9988 , 15.9994 , 39.948 , 4.0026 ,1.00797 / 


TX = 444.3807 + .02385 * T - 392.8292 * EXP ( -.0021357 * T ) 
A2 = 2. * (T-TX) / PI 
TX_T0 = TX - TO 
T1 = 1.9 *TX_T0/ 35. 

T3 = -1.7 *TX_T0/(35.**3) 

T4 = -0.8 * TX_T0 / ( 35. **4 ) 

TZ = TEMP ( Z , TX , T1 , T3 , T4 , A2 ) 


C** SECTION 1 
C** 

A = 90. 

D = AMIN 1 (Z, 105. ) 

C Integrate gM/T from 90 to minimum of Z or 105 km :- 
CALL GAUSS ( A, D, 1, R, TX , T1 , T3 , T4 , A2 ) 
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C The number 2. 1926E-8 = density x temperature/mean molecular weight at 90 km. 

EM = MOL_WT ( D ) 

TD = TEMP ( D , TX , T1 , T3 , T4 , A2 ) 

DENS = 2.1926E-8 * EM * EXP( -R / RGAS ) / TD 

FACTOR = AV * DENS 
PAR = FACTOR / EM 
FACTOR = FACTOR / 28.96 


C For altitudes below and at 105 km calculate the individual specie number 
C densities from the mean molecular weight and total density. 


IF ( Z. LE. 105) THEN 

DL = ALOG 1 0 ( DENS ) 

AN = ALOG 10 ( QN * FACTOR ) 

A A = ALOG 10 ( QA * FACTOR ) 

AHE = ALOG 10 ( QHE * FACTOR ) 

AO = ALOG 10 ( 2. * PAR * ( 1. -EM / 28.96 ) ) 

A02 = ALOG 10 ( PAR * ( EM * ( 1 .+Q02 ) / 28.96-1 . ) ) 
AH =0. 

C 

C** Return to calling program 
C 

RETURN 


END IF 


C** SECTION 2 : This section is only performed for altitudes above 105 km 

Q ** 


C Note that having reached this section means that D in section 1 is 105 km. 

C — 

C Calculate individual specie number densities from the total density and mean 
C molecular weight at 105 km altitude. 

DI(1) = QN * FACTOR 

DI(2) = PAR * (EM * ( 1 .+Q02) / 28.96- 1 .) 

DI(3) = 2. * PAR * (1 .- EM / 28.96) 

DI(4) = QA * FACTOR 
DI(5) = QHE * FACTOR 

C Integrate g/T from 105 km to Z km :- 
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CALL GAUSS ( D, Z, 2, R, TX , T1 , T3 , T4 , A2 ) 

DO 41 1=1,5 

DIT(I) = DI(I) * ( TD / TZ ) **( 1 ,+ALPHA(I)) * EXP( -EI(I) * 
& R / RGAS) 

IF ( DIT(I). LE. 0. ) DIT(I) = 1 .E-6 
41 CONTINUE 


C** This section calculates atomic hydrogen densities above 500 km altitude. 
C** Below this altitude , H densities are set to 10**-6. 

C** SECTION 3 
C** 

IF ( Z .GT. 500. ) THEN 

A1 =500. 

S = TEMP ( A 1 , TX , T1 , T3 , T4 , A2 ) 

DI(6) = 10.** (73.13 - 39.4 * ALOG10 (S) + 5.5 * ALOGIO(S) 

& *ALOG10(S)) 

CALL GAUSS ( Al, Z, 7, R, TX , T1 , T3 , T4 , A2 ) 

DIT(6) = DI(6) * (S/TZ) * EXP ( -EI(6) * R / RGAS ) 

ELSE 

DIT (6) = 1.0 

END IF 


C For altitudes greater than 105 km , calculate total density and mean 
C molecular weight from individual specie number densities. 

DENS=0 
DO 42 I = 1 , 6 

DENS = DENS + EI(I) * DIT(I) / AV 
42 CONTINUE 

EM = DENS * AV / ( DIT(1)+DIT(2)+DIT(3)+DIT(4)+DIT(5)+DIT(6) ) 
DL = ALOG10 (DENS) 

AN = ALOG10(DIT(1)) 
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A02 = ALOG10(DIT(2)) 

AO = ALOG10(DIT(3)) 

A A = ALOG10(DIT(4)) 

AHE = ALOG10(DIT(5)) 

AH = ALOG10(DIT(6)) 


RETURN 

END 


FUNCTION TEMP ( ALT , TX , T1 , T3 , T4 , A2 ) 

£*********** ****** ****** ************ He*****!!!********************************** 

He* 

C** Function subprogram 'TEMP' calculates the temperature at altitude ALT ** 

C** using equation (10) for altitudes between 90 and 125 km and equation ** 

C** (13) for altitudes greater than 125 km , from SAO Report 313. ** 

£** H=* 

C** Written by Mike Hickey, USRA ** 

Q* ********** ****** ********************************************** ************* 
** 


PARAMETER BB = 4.5E-6 


U = ALT - 125. 

IF ( U .GT . 0. ) THEN 

TEMP = TX + A2 * ATAN ( T1 * U * ( 1 . + BB * (U**2.5)) / A2 ) 
ELSE 

TEMP = TX + T 1 * U + T3 * (U**3) + T4 * (U**4) 

END IF 
END 


REALM FUNCTION MOL_WT ( A ) 

^* ************************************************ ********5^^^ C ****** ********** 
** 

£** ** 

C** Subroutine 'MOL_WT calculates the molecular weight for altitudes ** 

C** between 90 and 105 km according to equation (1) of SAO report 313. ** 


A24 



MARSHALL ENGINEERING THERMOSPHERE MODEL PROGRAM LISTING 

(PAGE 24 OF 28) 


C** Otherwise, MOL_WT is set to unity. ** 

£** ** 

C** Written by Mike Hickey, USRA ** 

£**************************************************************************** 

** 


DIMENSION B(7) 

DATA B / 28.15204 , -0.085586, 1.284E-4, -1.0056E-5, -1.021E-5, 
1 .5044E-6, 9.9826E-8 / 

IF ( A. GT. 105. ) THEN 


MOL_WT= 1. 


ELSE 


U = A - 100. 

MOL_WT = B ( 1 ) 

DO 1 1 = 2,7 

MOL_WT = MOL_WT + B (I) * U ** ( 1-1 ) 
1 CONTINUE 
END IF 
END 


SUBROUTINE GAUSS ( Z1 , Z2 , NMIN , R , TX , T1 , T3 , T4 , A2 ) 

** 

C** Subdivide total integration-altitude range into intervals suitable for ** 

C** applying Gaussian Quadrature , set the number of points for integration ** 

C** for each sub-interval , and then perform Gaussian Quadrature. ** 

C** Written by Mike Hickey, USRA, NASA/MSFC, ED44, July 1988. ** 

C************************ ********** ****************************************** 
** 

REAL*4 ALTMIN (9) , C(8,6), X(8,6), MOL_WT 
INTEGER NG (8) , NGAUSS , NMIN , J 

GRAVITY ( ALTITUDE ) = 9.80665 / ( ( 1 . + ALTITUDE / 

& 6.356766E3 )**2 ) 
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DATA ALTMIN / 90., 105., 125., 160., 200., 300., 500., 1500., 
& 2500. / 

DATA NG / 4 , 5, 6, 6, 6, 6, 6, 6/ 

C Coefficients for Gaussian Quadrature ... 


DATA C / .5555556 , .8888889 , .5555556 , .0000000 , ! n=3 
.0000000 , .0000000 , .0000000 , .0000000 , ! n=3 
.3478548 , .6521452 , .6521452 , .3478548 , ! n=4 
.0000000 , .0000000 , .0000000 , .0000000 , ! n=4 
.2369269 , .4786287 , .5688889 , .4786287 , ! n=5 
.2369269 , .0000000 , .0000000 , .0000000 , ! n=5 
.171 3245 , .36076 1 6 , .4679 1 39 , .4679 139,! n=6 
.36076 1 6 , . 1 7 1 3245 , .0000000 , .0000000 , ! n=6 
. 1 294850 , .2797054 , .38 1 830 1 , .4 1 79592 , ! n=7 
.3818301 , .2797054 , .1294850 , .0000000 , ! n=7 
.1012285,. 2223810,. 3137067,. 3626838, ! n=8 
.3626838 , .3137067 , .2223810 , .1012285/ ! n=8 


C Abscissas for Gaussian Quadrature ... 


DATA X / -.7745967 , .0000000 , .7745967 , .0000000 , ! n=3 
.0000000 , .0000000 , .0000000 , .0000000 , ! n=3 
-.861 1363 , -.3399810 , .3399810 , .861 1363 , ! n=4 
.0000000 , .0000000 , .0000000 , .0000000 , ! n=4 
-.9061798 , -.5384693 , .0000000 , .5384693 , ! n=5 
.9061798 , .0000000 , .0000000 , .0000000 , ! n=5 
-.9324695 ,-.66 12094, -.2386 192, .2386192 ,! n=6 
.6612094 , .9324695 , .0000000 , .0000000 , ! n=6 
-.9491079 , -.7415312 , -.4058452 , .0000000 , ! n=7 
.4058452 , .7415312 , .9491079 , .0000000 , ! n=7 
-.9602899 , -.7966665 , -.5255324 , -. 1 834346 , ! n=8 
.1834346, .5255324, .7966665, .9602899 /! n=8 


R = 0.0 


DO 2 K = NMIN , 8 

NGAUSS = NG (K) 

A = ALTMIN (K) 

D =AMIN1 (Z2, ALTMIN (K+l)) 

RR =0.0 

DEL = 0.5 * ( D - A ) 

J = NGAUSS -2 

DO 1 1=1, NGAUSS 

Z = DEL * ( X(I,J) + 1. ) + A 

RR = RR + C(I,J) * MOL_WT(Z) * GRAVITY(Z) / TEMP ( Z,TX,T1 ,T3, 
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& T4,A2 ) 

1 CONTINUE 

RR = DEL * RR 
R = R + RR 

IF ( D .EQ. Z2 ) RETURN 

2 CONTINUE 

RETURN 

END 


SUBROUTINE SLV ( DEN , ALT , XLAT , DAY ) 

** 

C** Subroutine 'SLV' computes the seasonal-latitudinal variation of density ** 
C** in the lower thermosphere in accordance with L. Jacchia, SAO 332, 1971. ** 
C** This affects the densities between 90 and 170 km. This subroutine need ** 
C** not be called for densities above 170 km, because no effect is observed. ** 
£** ** 

C** The variation should be computed after the calculation of density due to ** 
C** temperature variations and the density ( DEN ) must be in the form of a ** 
C** base 10 log. No adjustments are made to the temperature or constituent ** 
C** number densities in the region affected by this variation. ** 


Q** 


** 


c** 

DEN = density (log 10) 


** 

C** 

ALT = altitude (km) 


** 

c** 

XLAT = latitude (rad) 


** 

c** 

C** 

DAY = day number 

** 

** 


£********************************* ******************************************* 
** 


C** initialize density (DEN) = 0.0 
C 

DEN = 0.0 


C 

C** check if altitude exceeds 170 km 
C 

IF ( ALT. GT. 170. ) RETURN 
C 

C** compute density change in lower thermosphere 
C 


Z = ALT - 90. 
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X = -0.0013 *Z*Z 
Y = 0.0172 * DAY + 1.72 
P = SIN (Y) 

SP = ( SIN (XLAT) ) **2 
S = 0.014 * Z * EXP (X) 

D = S * P * SP 
C 

C** check to compute absolute value of 'XLAT' 
C 

IF ( XLAT. LT. 0. ) D = -D 
DEN = D 


RETURN 

END 


SUBROUTINE SLVH ( DEN , DENHE , XLAT , SDA ) 

** 

C** Subroutine 'SLVH' computes the seasonal-latitudinal varaition of the ** 
C** helium number density according to L. Jacchia, SAO 332, 1971. This ** 
C** correction is not important below about 500 km. ** 

£** ** 

C** DEN = density (log 10) ** 

C** DENHE = helium number density (log 10) ** 

C** XLAT = latitude (rad) ** 

C** SDA = solar declination angle (rad) ** 

** 


DO = 10. ** DENHE 
A = ABS ( 0.65 * ( SDA / 0.40909079 ) ) 

B = 0.5 * XLAT 
C 

C** Check to compute absolute value of 'B' 

C 

IF ( SDA. LT. 0. ) B = -B 
C 

C** compute X, Y, DHE and DENHE 
C 

X = 0.7854 - B 
Y = ( SIN (X) ) ** 3 
DHE= A * ( Y - 0.35356 ) 

DENHE = DENHE + DHE 
C 
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C** compute helium number density change 
C 

D1 = 10. ** DENHE 
DEL= D1 - DO 
RHO= 10. ** DEN 
DRHO = ( 6.646E-24 ) * DEL 
RHO = RHO + DRHO 
DEN = ALOGIO (RHO) 

RETURN 

END 


SUBROUTINE FAIR5 ( DHEL1 ,DHEL2 ,DLG1 ,DLG2 ,IH .FDHEL .FDLG ) 


£****************************************************************************** 
C** This subroutine fairs between the region above 500 km, which invokes the ** 

C** seasonal-latitudinal variation of the helium number density ( subroutine ** 

C** SLVH ), and the region below, which does not invoke any seasonal- ** 

C** latitudinal variation at all. ** 

£** ** 

C** INPUTS: DHEL1 = helium number density before invoking SLVH ** 


c ** 

C** 

c** 

c** 

C** 


DHEL2 = helium number density after invoking SLVH ** 

DLG1 = total density before invoking SLVH ** 

DLG2 = total density after invoking SLVH ** 

IH = height ( km ) - INTEGER ** 

IBFH = base fairing height ( km ) — INTEGER ** 

C** OUTPUTS: FDHEL = faired helium number density ** 

C** FDLG = faired total density ** 

£** ** 

C** Written by Bill Jeffries, CSC, Huntsville, AL. ** 

C** ph. (205) 830-1000, x3 11 ** 

£****************************************************************************** 


DIMENSION CZ ( 6 ) 

PARAMETER IBFH = 440 

DATA CZ / 1.0, 0.9045085, 0.6545085, 0.3454915, 0.0954915, 0.0/ 
C Height index 

I = ( IH - IBFH ) / 1 0 + 1 
C Non-SLVH fairing coefficient 
CZI = CZ ( I ) 

C SLVH fairing coefficient 
SZI= 1.0-CZI 
C Faired density 

FDLG = ( DLG1 * CZI ) + ( DLG2 * SZI ) 

C Faired helium number density 

FDHEL = ( DHEL1 * CZI ) + ( DHEL2 * SZI ) 

RETURN 

END 
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USER'S SOFTWARE IMPLEMENTATION GUIDE 

AND 

PROGRAM LISTINGS 

A SIMULATION OF SMALL-SCALE 
THERMOSPHERIC DENSITY 
VARIATIONS FOR ENGINEERING 
APPLICATIONS 

(For orbital inclinations < 40°) 



USER'S SOFTWARE IMPLEMENTATION GUIDE 
FOR ORBITAL INCLINATIONS LESS THAN OR EQUAL TO 40° 

This section describes how to use the wave simulation 
model. In this section we employ the following conventions: 
subroutine and common blocks are written using all upper 
case letters, while only the first character of a variable 
name is written using upper case letters. There are several 
subroutines employed in the model. 'fheir names, order of 
execution and functions are provided in Table Bl, while 
their input and corresponding output are provded in Table 
B2 . They are also included in Program Listing on p. B9. 


Table Bl. Subroutines, Order of Execution and Functions 


Subroutine Name 

Input Parameters 

Output Parameters 

5 ETWAVES c 

Pickbins: user input 
F^bin and Apbin (if 
Pickbins true): user input 
IseedO: user input 
Warnings: user input 

IseedO: through argument 
list 

Fi3bin, Apbin and 
Warnings: all through 
COMMON /WAVE/ 

GETW a V ES 

Metden, Pi3in, apin, and 
IseedO: ail through 
argument list 

Totden: through argument 
list 

WAVES 

Apin, Fiyn, Fixdstep, Dt, 
Tmax, Iseed, Reset, 
Pickbins, F^bin, Apbin, 
Warnings: all through 
argument list 

Wave: through argument 
list 

GAUSSD 

N and Iseed: through 
argument list 

Noise: through argument 
list 

RANI 

Idtun (seed): through 
argument list 

Rani: function 
subprogram value 
returned 


These subroutines have been designed to be used in 
conjunction with both the MET model and an orbit generator, 
although in principle they will work with anu empirical 
thermospheric density model. The subroutine WAVES is called 
and controlled by the subroutine GETWAVES . The input and 
output of WAVES are thoroughly described in the comments 
accompanying that subroutine, and will not need further 
explanation here. The subroutine GAUSSD and the function 
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subprogram RANI, that together generate an array of Gaussian 
distributed random numbers, need no further explanation here 
except that the seed must be a large odd integer. The ref- 
erences from which these portable and reliable codes were 
obtained is given in the comments section of each of their 
program codes . 


Table B2 . Subroutines and Their Input and Output. 


Subroutine 

Name 

Execution Order 
and Call Type 

Function 

SETWAVES 

1 Single call 

Performs routine setup. Requires user input 
for option to explicitly select F 13 bins, the 
initial random number seed, and the option 
to suppress warning messages. 

GETWAVES 

2 Multiple call 

Handles bookkeeping for WAVES 
subroutine. Inputs through argument list. 
Calls WAVES subroutine. 

WAVES 

3 Multiple call 

Computes array of stochastic wave 
amplitudes. Inputs through argument list 
Calls GAUSSD subroutine. 

GAUSSD 

4 Multiple call 

Computes array of Gaussian random 
numbers. Inputs through argument list 
rails RANI subroutine 

RANI 

5 Multiple call 

Computes a single uniform random deviate 
between 0 and 1 . Inputs through argument 
list 


Correct implementation of the set of subroutines will 
require some additional comments regarding the subroutines 
SETWAVES and GETWAVES . These are now briefly described. 


Subroutine SETWAVES 

The subroutine SETWAVES must be called initially to set 
some of the inputs. If the decisions made during the execu- 
tion of SETWAVES are adhered to for the entire simulation or 
sets of simulations , SETWAVES need not be called again. In 
this case Iseed) , Pickbins and Warnings, as well as F^bin 
and Apbin if Pickbins is true, will all remain unchanged. 

If, however, Pickbins or Warnings need to be changed, then 
SETWAVES will need to be recalled. In practice it is un- 
likely, for example, that any user would want to run a simu- 
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lation first without suppressing warning statements and then 
later suppressing them. It would seem equally unlikely that 
during the course of a set of simulations any user would 
need to change from using a predetermined set of ap and F ^3 
values to using predetermined a p and F 13 bin values. 

Subroutine GETWAVES 

The subroutine GETWAVES needs to be called directly 
after the total neutral mass density, that has been 
calculated from the mean-state thermospheric model 
(e . g . , MET ) , has been returned to the orbit generator 
program. GETWAVES must be called every time that a new 
density value has been returned from the mean-state density 
model . 

The first input parameter in the arument list of 
GETWAVES is Metden, the mean-state total neutral mass 
density output from the standard MET mode. An output Metden 
is modified by the addition of a stochastic component, as 
described earlier in the report. The other parameters in 
the argument list are F^in, A p in and IseedO. As described 
previously, if SETWAVES is called once only, IseedO will 
remain unchanged. the input parameters F 33 in and A p in must 
be available within the orbit generator program whenever 
Pickbins is false so that they can be properly input to 
GETWAVES (in that case WAVES calculates the bins internally 
using the F 13 in and A p in values) . If Pickbins is true, so 
that the bins are explicitly user-selected, F 13 in and A D in 
are not used at all. p 

Other parameters are passed to GETWAVES through a Common 
block (WAVE). The Common block must be included within the 
orbit generatore program exactly as specified in GETWAVES. 

The parameters in this Common block include the four logical 
constants Fixdstep, Reset, Pickbins and Warnings. The last 
three of these have their values set in SETWAVES (as 
previously discussed) . The first logical constant, 

Fixdstep, must be set within the orbit generator program. 

Its setup belongs there, and not within SETWAVES, because 
the orbit generator program usually handles the time step, 
dt • For a single simulation, Fixdstep must remain constant. 
It can be changed only when Reset is true. 

The remaining parameters in the Common block are Dt, Tmax, 
Iseed, F 33 bin and A p bin. As already discussed, F 33 bin and 
Apbin may have been previously set in SETWAVES. The time 
step, Dt, and the simulation time, Tmax, are set within the 
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orbit generator program. If Reset is true during the course 
of a simulation, the random number seed, Iseed, is reset by 
adding the integer 2 to its previous value. This method 
will ensure that Iseed remains large and odd while not 
interrrupting program execution (as would occur is a new 
value of IseedO was required to be manually input from a 
terminal by the user) > 

A single simulation should not be allowed to exceed 
Tmax. If the elapsed time exceeds Tmax, Reset is reset to 
true, Elaptime is reset to zero, the random number seed, 
Iseed, is updated and a new simulation (which may follow the 
previous one) is initiated. Similarly, a new simulation is 
initiated if Pickbins is false and the ap value changes. We 
don't expect F]_ 3 to change during the course of these 
relatively short duration simulations (less than a day), so 
no option exists for that remote possibility. If Tmax 
exceeds three hours and Pickbins is false, then it is likely 
that ap will change because it is a three-hour index. 
Therefore, the simulations are not allowed to exceed three 
hours. If the user inputs a value of Tmax that is greater 
than three hours, the subroutine WAVES will reset that value 
to three hours and issue a warning (if Warnings is true) . 

Modifications Required In User-Supplied Orbit Generator 
Progam 

The five following parameters need to be set by the 
user directly before calling GETWAVES within the user- 
written orbit generator program: 


a) 

Dt 

(set 
step ) 

once 

only if using a fixed time 

b) 

Fixdstep 

(set 

once 

only) 

c) 

F i3 in 

( set 

once 

only) 

d) 

Apin 

( set 

once 

only) 

e) 

Tmax 

( set 

once 

only) _ 


The subroutines that need to be called from within the 
user-written orbit generator program and their calling 
sequence are: 

a) Thermospheric density model (e.g. , MET) 

then immediately call next routine once only. 

b) SETWAVES 

then immediately call next routine multiple 
times 

c) GETWAVES 
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PROGRAM LISTINGS 


The following pages contain the program listings for 
the subroutines SETWAVES, GETWAVES and GAUSSD and the 
function subprogram RANI. The coding for this software is 
written in FORTRAN 77. It is not machine specific and 
should, therefore, be portable between different machines. 
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SUBROUTINE SETWAVES (ISEEDO) 


IMPLICIT NONE 

INTEGER APBIN, F13BIN, I SEED, ISEEDO 

REAL DT, TMAX 

CHARACTER* 1 ANSWER 

LOGICAL PICKBINS, WARNINGS, FIXDSTEP, RESET 

COMMON / WAVE / FIXDSTEP, DT, TMAX, ISEED, RESET, PICKBINS, 

F13BIN, APBIN, WARNINGS 

Z Set option for explicitly selecting low or high F10.7 and ap bins for 
Z wave simulator. The two F10.7 bins relate to the 13-month smoothed 
2 value of F10.7 lying below 138 (low bin) or above 138 (high bin). The 
z two ap bins correspond to ap values below the 50th percentile (the 
Z low ap bin for ap values less than or equal to 7) and those above the 
: 50th percentile (the high ap bin for ap greater than 7) . If the option 

j of explicitly selecting these bins is rejected, the bins are chosen 
C according to the values of 3 -hourly ap and the 162-day mean F10.7 
; (assumed to approximate the ’13 -month smoothed FI 0.7) that are input 
2 to the MET model. The option is particularly important if 13 -month 
Z smoothed values of predicted ap values are employed in simulations, 

: as they generally have values that lay solely in the **** ap bin. 

PRINT *,' Explicitly select F10 and ap bins for wave simulator?' 
PRINT *, ' (Y/N) ' 

1 FORMAT (Al) 

READ 1, ANSWER 

IF (ANSWER. EQ. 'Y' .OR. ANSWER. EQ. 'y') THEN 
PICKBINS = .TRUE. 

PRINT * , 'Low (=1) or high (=2) 13-month F10.7 bin?' 

READ *, F13BIN 

PRINT * , ' Low (=1) or high (=2) ap bin?' 

READ *, APBIN 
ELSE 

PICKBINS » .FALSE. 

END IF 

0 An initial seed for the random number generator used in the wave 
C simulator is required. For subsequent simulations the seed is 
: internally incremented by 2. 

PRINT *,'Set initial wave simulator seed (large odd integer)' 
READ *, ISEEDO 

: Warnings will be typed from the wave simulator if certain parameters 

C acre incorrectly specified. For am "expert" user these cam be 
suppressed by setting the warnings option to false. 

PRINT *, 'Suppress warning or error messages? (Y/N)' 

READ 1, ANSWER 

IF (ANSWER. EQ. 'Y' .OR. ANSWER. EQ. 'y') THEN 
WARNINGS - .FALSE. 

ELSE 

WARNINGS - .TRUE. 

END IF 

C End setup for wave simulator 


RETURN 

END 
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SUBROUTINE 

IMPLICIT 

INTEGER 

REAL 

• 

LOGICAL 


GETWAVES (METDEN, FI 3 IN, APIN, ISEEDO, TOTDEN) 

N ONE 

JSTE P, NRESET, I SEED, ISEEDO, IWAVE, F13BIN, APBIN 
MET DEN, FI 3 IN, APIN, TOTDEN, DT, TIMESTEP (720) , 
ELAP TIME , TMAX, OLDAP, WAVE (720) 

FIXDSTEP, RESET, PICKBINS, WARNINGS 


COMMON / WAVE / FIXDSTEP, DT, TMAX, ISEED, RESET, PICKBINS, 

F13BIN, APBIN, WARNINGS 


IF (ELAPTIME. GT.TMAX .OR. (.NOT. RESET .AND. .NOT. PICKBINS 
. .AND. OLDAP. NE. APIN) ) RESET - .TRUE. 

IF (RESET) THEN 
NR ESE T - NRESET+1 

ISEED - ISEEDO+2 * NRESET ! this keeps seed odd 

OLDAP « APIN 
ELAPTIME » 0.0 
END IF 


Note that WAVES is called once only for fixed time steps 
IF (RESET. OR. .NOT. FIXDSTEP) THEN 
IWAVE - 0 

CALL WAVES (APIN, F13IN, FIXDSTEP, DT, TMAX, ISEED, RESET, 
PICKBINS, F13BIN, APBIN, WARNINGS, WAVE) 

END IF 

IWAVE « IWAVE+1 

ELAPTIME « ELAPTIME+DT 

TOTDEN * ( 1 . 0+WAVE ( IWAVE) ) * METDEN 

RETURN 

END 
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SUBROUTINE WAVES (APIN, F13IN, FIXDSTEP, DT, TMAX, ISEED, RESET, 

PICKBINS, F13BIN, APBIN, WARNINGS, WAVE) 

Z This subroutine /model is based on the analysis of AEE NACE data (see 
z accompanying report by Michael P. Hickey) . 

2 it is valid for low- inclination, approximately circular orbits 

3 between about 250 and 500 km altitude. The time resolution of the 

z original data and for this model is 15.0 seconds. The perturbations 
Z output by this subroutine are used to modify the total mass density 
Z of the MET model and are intended for simulations related to 
Z atmospheric drag effects (e.g., guidance, navigation and control). 

Z INPUTS: APIN is the same ap value that is input to MET model; 

F13IN is the 13 -month smooth value of monthly mean F10.7; 

FIXDSTEP tells subroutine if time steps will be constant or 
if they will vary during the course of simulation; 

DT is either the fixed time step if FIXDSTEP is true 

or is the next time step if FIXDSTEP is false. In 
either case , the minimum allowable value of DT is 
15 second^ (the time resolution of original data) . 

If the input value of DT is less than 15 s it is 
automatically reset to 15.0; 

TMAX is the maximum time duration (sec) of simulation. 

It is advised that TMAX not exceed 3 hours (the 
time resolution of the ap index) ; 

I SEED is a large odd integer for random number generator; 

RESET tells subroutine to go back to setup mode if it is 

true. It should be set false on every subsequent 
call to WAVES within a single simulation. It should 
be set to true every three hours of simulation or 
whenever ap or F13 vary; 

PICKBINS tells subroutine to override APIN and F13 IN inputs 
and explicitly select required AP and F13 bins ; 
F13BIN is the explicitly selected F13 bin when PICKBINS is 

true; 

APBIN is the explicitly selected Ap bin when PICKBINS is 

true; 

WARNINGS tells subroutine to type all warnings that occur 
when either invalid input parameters sure reset 
internally or when WAVES should be recalled in 
setup mode. If WARNINGS is false, no warnings are 
typed to the screen , although the user should be 
aware that internal resetting may be occuring if 
the subroutine is being misused. 

OUTPUT: WAVE is the wave percentage amplitude expressed as a 

decimal. WAVE is an array dependant on FIXDSTEP. 

For FIXDSTEP false, only the first element of WAVE 
is nonzero; subroutine WAVES is then called for 
every time step to determine WAVE(l) . 

For FIXDSTEP true, WAVE has l+TMAX/15.0 elements 
for DT less than or equal to 15s and 1+TMAX/DT 
elements for DT greater than 15s. 

v. Subroutine WAVES is then called once only to 

C determine WAVE (i) , i-1, l+min(DT,15.0) . 

Additional Subroutines Employed: GAUSSD and RANI (random / generators) 
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IMPLICIT 

INTEGER 

REAL 


LOGICAL 

• 

DATA 

DATA 

DATA 

DATA 


NONE 

ISEED, F13BIN, APBIN, N, K, I 
TMAX, 

APIN, F13IN, DT, WAVE(1+TMAX/DT) , NOISE (721) , PHI, 
STNDEV, AMP (721) , RATIO, ELTIME, PHI1, PHI2 , 

MARK OV ( 2 ) , SCALE (2,2) , AR2(4) 

FIXDSTEP, SETUP, RESET, PICKBINS, WARNINGS, 
STOPFLAG 

MARKOV / 0.8531, 0.8324 / 

AR2 / 0.4672, 0.2841, 0.4902, 0.3046 / 

S CALE / 2.00E-2, 1.80E-2, 1.91E-2, 1.92E-2 / 

SETUP / .TRUE. / 


IF (RESET) SETUP * .TRUE. 


Check that DT and TMAX are both nonzero. Type fatal error message 
and terminate program execution if either one of them is zero. 
STOPFLAG * .FALSE. 

IF (DT.EQ.O.) THEN 
PRINT 1 

STOPFLAG * .TRUE. 

ELSE IF (TMAX.EQ.O.) THEN 
PRINT 2 

STOPFLAG * .TRUE. 

END IF 

IF (STOPFLAG) STOP 


Further check that DT is no smaller than 15 sec and that TMAX is no 
greater than 3 hrs. If they are, reset them to their nominal values. 
IF (DT.LT.15.0) THEN 
DT - 15.0 

IF (WARNINGS) PRINT 3 
END IF 

IF (TMAX.GT. 1. 08 E4) THEN 
TMAX - 1.08E4 
IF (WARNINGS) PRINT 4 
END IF 


Commence setting F13 & ap bins, white noise array and wave amplitude 
array on a constant 15-sec time grid. 

IF (SETUP) THEN 

SETUP - .FALSE. 

RESET » .FALSE. 

If PICKBINS is true, then ap and F13 bins are directly input. 
Otherwise, they sure calculated from the input values of ap and F13. 
IF ( .NOT. PICKBINS) THEN 
IF (APIN.LE. 7.0) THEN 
APBIN « 1 
ELSE 

APBIN * 2 
END IF 

IF (F13IN.LE. 138.0) THEN 
F13BIN - 1 
ELSE 

F13BIN - 2 
END IF 
END IF 


BIO 



Set phi and estimated white noise standard deviation of process 
IF (F13BIN.EQ.1) THEN 
PHI - MARKOV (APBIN) 

ELSE 

PHI1 - AR2 (2*APBIN-1) 

PHI 2 « AR2 (2*APBIN) 

END IF 

STNDEV * SCALE (F13BIN, APBIN) 

N * l+TMAX/15.0 

Generate white noise of zero mean and unit standard deviation 
CALL GAUSSD (NOISE, N, I SEED) 

Generate array of wave amplitudes with 15-sec time resolution after 
first initializing array to zero. 

DO I * 1, N 
AMP (I) = 0.0 
END DO 

AMP ( 1) * NOISE (1) 

IF (F13BIN.EQ.1) THEN 
DO I * 2, N 

AMP (I) = PHI*AMP (I — 1 ) + NOISE (I) 

END DO 
ELSE 

AMP ( 2 ) = PHI1*AMP ( 1) + NOISE (2) 

DO I = 3, N 

AMP (I) = PHI1*AMP (1-1) + PHI2*AMP (I— 2 ) + NOISE (I) 

END DO 
END IF 

Now make variance fit data by increasing standard deviation from 
unity to STNDEV. 

DO I * 1, N 

AMP (I) = STNDEV* AMP ( I ) 

END DO 

If FIXDSTEP is true transform from 15-s time resolution array to DT 
time resolution array. . . 

IF (FIXDSTEP) THEN 
RATIO * DT/15.0 
DO 1*1, 1+TMAX/DT 

K * 1+INT (RATIO*REAL(I-l) ) 

WAVE (I) * AMP (K) 

END DO 

. . . else keep track of elapsed time 
ELSE 

ELTIME » DT 

K - 1+INT (ELTIME/ 15.0) 

WAVE ( 1 ) * AMP (K) 

END IF 
RETURN 

Next section is not setup, and should only be executed for FIXDSTEP 
false . 

ELSE 

IF ( .NOT. FIXDSTEP) THEN 
ELTIME - ELTIME+DT 
IF (ELTIME. GT.TMAX) THEN 
PRINT 5 
WAVE ( 1 ) =» 


0.0 
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ELSE 

K - 1+INT ( ELTIME/ 15.0) 

WAVE ( 1 ) - AMP (K) 

END IF 
RETURN 
ELSE 

IF (WARNINGS) PRINT 6 
END IF 
END IF 

FORMAT ( ' FATAL ERROR: DT is zero in WAVES') 

FORMAT (' FATAL ERROR: TMAX is zero in WAVES') 

FORMAT (' WARNING: DT has been increased to 15 seconds') 

FORMAT (' WARNING: TMAX has been decreased to 3 hours') 

FORMAT (' WARNING: Elapsed time > TMAX. WAVE now set to 0.0' , /, 

• Recall WAVES in setup mode') 

FORMAT (' WARNING: Fixed time step. Recall WAVES in setup mode') 

RETURN 

END 
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SUBROUTINE GAUSSD (X, N, ix) 

routine to produce vector of gaussian distributed random nu mber s 
in X(n) . The series in X(n) will have approximate zero mean and 
standard deviation of 1.0. Polar Method algorithm taken from 
p. 104, Seminumerical Algorithms, 

vol 2 of The Art of Computer Programming by D. E. Knuth, 1969. 

ix is random number seed , i . e . , large odd integer . 

implicit none 
integer i, ix, iy, n 
REAL X(N) , RANI 
real s, vl, v2 


I - 0 
iy = ix 

do while (i .It. n) 

Vl = 2.0*ranl(iy) - l.o 
V2 = 2 . 0*ranl (iy) - 1.0 

S = V1*V1 + V2*V2 

if (s .It. 1.0) then 

S - SQRT (—2.0 *ALOG ( S ) / S) 

1 = 1 + 1 
X(I) = Vl * s 
1 = 1 + 1 

if (i .le. n) then 
X(I) = V2 * S 
end if 

end if 

end do 

ix = iy 
RETURN 

END 
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FUNCTION RANI (IDUM) 

Returns a uniform random deviate between 0.0 and 1.0. Set IDUM to any 
negative value to initialize or reinitialize the sequence. 

Copied from Numerical Recipes , W.H. Press, B.P. Flannery, 

S.A. Teukolsky & W.T. Vetterling, Cambridge Univ. Press, 1986, p.196 
IMPLICIT NONE 

INTEGER IFF, 1X1, 1X2, 1X3, J, Ml, M2, M3, IA1, IA2 , IA3 , 

IC1, IC2, IC3, IDUM 
REAL RM1, RM2 , R, RANI 

DIMENSION R(97) 

PARAMETER (Ml-259200, IA1-7141, IC1-54773, RM1-1./M1) 

PARAMETER (M2-134456, IA2-8121, IC2-28411, RM2-1./M2) 

PARAMETER (M3-243000, IA3-4561, IC3-51349) 

DATA IFF / 0/ 1 as above, initialize on first call even if 

I IDUM is not negative 
IF (IDUM.LT.O .OR. IFF.EQ.O) THEN 
IFF-1 

1X1— MOD ( IC1— IDUM , Ml ) l Seed the first routine, 

IX1-MOD ( IA1*IX1+IC1 , Ml ) 

1X2— MOD ( 1X1 , M2 ) ! emd use it to seed the second 

IX1-MOD ( IA1*IX1+IC1 , Ml ) 

1X3— MOD ( 1X1 , M3 ) ! and third routines. 

DO J - 1, 97 ! Fill the table with sequential 

1X1— MOD (IA1*IX1+IC1 ,M1) ! uniform deviates generated by 

1X2— MOD (IA2 *1X2+1 C2 ,M2 ) 1 the first two routines 

R ( J) = ( FLOAT ( 1X1 ) + FLOAT ( 1X2 ) *RM2 ) *RM1 ! Low and high- 

END DO ! order pieces combined here. 

IDUM—1 
END IF 

1X1— MOD ( IA1*IX1+IC1 , Ml ) ! Except when initializing, this is 
1X2— MOD(IA2*IX2+IC2,M2) ! where we start. Generate the next 
1X3— MOD ( IA3*IX3+IC3 , M3 ) ! number for each sequence. 

J— 1+ (97*1X3 ) /M3 1 Use the third sequence to get an 

! integer between 1 and 97 
IF ( J. GT . 97 .OR. J.LT.l) PAUSE 

RANl-R(J) I return that table entry, 

R ( J) -( FLOAT (IX1)+FL0AT( 1X2 )*RM2)*RM1 I and refill it. 

RETURN 

END 
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APPENDIX C 

USER'S SOFTWARE IMPLEMENTATION GUIDE 

AND 

PROGRAM LISTINGS 


AN ENGINEERING MODEL FOR THE SIMULATION OF SMALL-SCALE 
THERMOSPHERIC DENSITY VARIATIONS FOR ORBITAL INCLINATIONS 

GREATER THAN 40° 



USER’S SOFTWARE IMPLEMENTATION GUIDE 
FOR ORBITAL INCLINATIONS GREATER THAN 40° 


This section describes how to use the wave simulation model. In this section we em- 
ploy the following conventions: subroutines and common blocks are written using all up- 
percase letters, while only the first character of a variable name is written using uppercase 
letters. There are several subroutines employed in the model. Their names, order of exe- 
cution and functions are provided in Table Cl, while their input and corresponding output 
are provided in Table C2. 


Table Cl. Subroutines, order of execution and functions. 


Subroutine 

Name 

Execution Or- 
der & call type 

Function 

SETWAVE2 

1 Single call 

Performs routine setup. Requires user in- 
put for option to explicitly select ap bin, 
the initial random number seed, and the 
option to suppress warning messages. 

GETWAVE2 

2 Multiple call 

Handles bookkeeping for WAVES sub- 
routine. Inputs through argument list. 
Calls WAVES subroutine. 

WAVES2 

3 Multiple call 

Computes stochastic wave amplitude. In- 
puts through argument list. 

Calls GAUSSD subroutine. 

GAUSSD 

4 Multiple call 

Computes array of Gaussian random 
numbers. Inputs through argument list. 
Calls RANI subroutine 

RANI 

5 Multiple call 

Computes a single uniform random devi- 
ate between 0 and 1 . Inputs through ar- 
gument list. 


These subroutines have been designed to be used in conjunction with both the MET model 
and an orbit generator, although in principle they will work with any empirical 
thermospheric density model. The subroutine WAVES2 is called and controlled by the 
subroutine GETWAVE2. The input and output of WAVES2 are thoroughly described in 
the comments accompanying that subroutine, and will not need further explanation here. 
The subroutine GAUSSD and the function subprogram RANI, that together generate an 
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array of Gaussian distributed random numbers, need no further explanation here except 
that the seed must be a large odd integer. The references from which these portable and 
reliable codes were obtained is given in the comments section of each of their program 
codes. 


Table C2. Subroutines and their input and output. 


Subroutine Name 

Input Parameters 

Output Parameters 

SETWAVE2 

Pickbins: user input 
Apbin (if Pickbins true): 
user input 
IseedO: user input 
Warnings: user input 

IseedO: through argument 
list 

Apbin and Warnings: all 
through COMMON 
/WAVE/ 

GETWAVE2 

Metden, latin, apin, and 
IseedO: all through argu- 
ment list 

Totden: through argument 
list 

WAVES2 

Apin, latin, Dt, Tmax, 
Iseed, Reset, Pickbins, 
Apbin, Warnings: all 
through argument list 

Wave: through argument 
list 

GAUSSD 

N and Iseed: through ar- 
gument list 

Noise: through argument 
list 

RANI 

Idum (seed): through ar- 
gument list 

Rani: function subprogram 
value returned 


Correct implementation of the set of subroutines will require some additional com- 
ments regarding the subroutines SETWAVE2 and GETWAVE2. These are now briefly 
described. 

Subroutine SETWAVE2 

The subroutine SETWAVE2 must be called initially to set some of the inputs. If the 
decisions made during execution of SETWAVE2 are adhered to for the entire simulation 
or sets of simulations, SETWAVE2 need not be called again. In this case IseedO, Pickbins 
and Warnings, as well as apbin if Pickbins is true, will all remain unchanged. If, however, 
Pickbins or Warnings need to be changed, then SETWAVE2 will need to be recalled. In 
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practice it is unlikely, for example, that any user would want to run a simulation first with- 
out suppressing warning statements and then later suppressing them. It would seem 
equally unlikely that during the course of a set of simulations any user would need to 
change from using a predetermined a p value to using a predetermined a p bin values. 

Subroutine GETWAVE2 

The subroutine GETWAVE2 needs to be called directly after the total neutral mass 
density, that has been calculated from the mean-state thermospheric model (e g., MET), 
has been returned to the orbit generator program. GETWAVE2 must be called every time 
that a new density value has been returned from the mean-state density model. 

The first input parameter in the argument list of GETWAVE2 is Metden, the mean- 
state total neutral mass density output from the standard MET model. On output Metden 
is modified by the addition of a stochastic component, as described earlier in this report. 
The other parameters in the argument list are latin, a p in and IseedO. As described previ- 
ously, if SETWAVE2 is called once only, IseedO will remain unchanged. The input pa- 
rameters latin and a p in must be available within the orbit generator program whenever 
Pickbins is false so that they can be properly input to GETWAVE2 (in that case WAVES2 
calculates the bins internally using the latin and a p in values). If Pickbins is true, so that the 
a p bin is explicitly user-selected, a p in is not used at all. We strongly emphasize that a lati- 
tude bin cannot be explicitly selected even if Pickbins is true, because latitude will always 
vary (and the latitude bins will too) around an orbit for which this simulator was designed 
(i.e., for inclinations grater than about 45°). 

Other parameters are passed to GETWAVE2 through a Common block (WAVE). 
This Common block must be included within the orbit generator program exactly as 
specified in GETWAVE2. The parameters in this Common block include the three logical 
constants Reset, Pickbins and Warnings. They have their values set in SETWAVE2 (as 
previously discussed). 

The remaining parameters in the Common block are Dt, Tmax, Iseed and a p bin. As 
already discussed, a p bin may have been previously set in SETWAVE2. The time step, Dt, 
and the simulation time, Tmax, are set within the orbit generator program. If Reset is reset 
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to true during the course of a simulation, the random number seed, Iseed, is reset by add- 
ing the integer 2 to its previous value. This method will ensure that Iseed remains large 
and odd while not interrupting program execution (as would occur if a new value of 
IseedO was required to be manually input from a terminal by the user). 

A single simulation should not be allowed to exceed Tmax. If the elapsed time exceeds 
Tmax, Reset is reset to true, Elaptime is reset to zero, the random number seed (Iseed) is 
updated, and a new simulation (which may follow the previous one) is initiated. Similarly, 
a new simulation is initiated if Pickbins is false and the a p value changes. If Tmax exceeds 
three hours and Pickbins is false, then it is likely that a p will change because it is a three 
hour index. Therefore the simulations are not allowed to exceed three hours. If the user 
inputs a value of Tmax that is greater than three hours, the subroutine WAVES2 will reset 
that value to three hours and issue a warning (if Warnings is true). 

Modifications Required in User-Supplied Orbit Generator Program 

The four following parameters need to be set by the user directly before calling GET- 
WAVES within the user-written orbit generator program: 


a) Dt 

(set once only if using a fixed time step) 

b) latin 

(set on every call) 

c) a p in 

(set once only) 

d) Tmax 

(set once only) 


The subroutines that need to be called from within the user-written orbit generator pro- 
gram and their calling sequence are 

a) Thermospheric density model (e g., MET) 

then immediately call next routine once only 

b) SE1WAVE2 

then immediately call next routine multiple times 

c) GETWAVE2 
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PROGRAM LISTINGS 


The following pages contain the program listings for the subroutines SETWAVE2, 
GETWAVE2, WAVES2 and GAUSSD and the function subprogram RANI. The coding 
for this software is written in FORTRAN 77. It is not machine specific and should there- 
fore be portable between different machines. 



n n n n n n 


SUBROUTINE SETWAVE2 (ISEEDO) 


IMPLICIT 

INTEGER 

REAL 

CHARACTER* 1 
LOGICAL 


NONE 

APB IN, LATBIN, I SEED, ISEEDO 

DT, TMAX 

ANSWER 

PICKBINS, WARNINGS, RESET 


COMMON / WAVE / DT, TMAX, I SEED, RESET, PICKBINS, 

> APB IN, WARNINGS 

C...Set option for explicitly selecting low or high ap bin for wave 
C. . .simulator . 

Z. . .Two ap bins correspond to ap values below the 50th percentile (the 
Z . . .low aD bin for ap values less than or equal to 7) and those above the 
I... 50th percentile '(the high ap bin for ap greater than 7) . If the option 
Z. . .of explicitly selecting these bins is rejected, the bin is chosen 
Z. . .according to the values of 3 -hourly ap that are input to the MET model. 
PRINT Explicitly select ap bin for wave simulator?' 

PRINT *, ' (Y/N) ' 

1 FORMAT (Al) 

READ 1, ANSWER 

IF (ANSWER. EQ. 'Y' .OR. ANSWER. EQ . 'y' ) THEN 
PICKBINS = .TRUE. 

PRINT * , ' Low (=1) or high (=2) ap bin?' 

READ *, APB IN 
ELSE 

PICKBINS = .FALSE. 

END IF 


. . .An initial seed for the random number generator used in the wave 
...simulator is required. For subsequent simulations the seen is 
. . . internally incremented by 2 . 

PRINT *,'Set initial wave simulator seed (large odd integer)' 
READ * , ISEEDO 

. . .Warnings will be typed from the wave simulator if certain parameters 
. . .are incorrectly specified. For an "expert" user these can be 
...suppressed by setting the warnings option to false. 

"PRINT *, 'Suppress warning or error messages? (Y/N)' 

READ 1, ANSWER 

IF (ANSWER. EQ. 'Y' .OR. ANSWER. EQ. 'y' ) THEN 
WARNINGS = .FALSE. 

ELSE 

WARNINGS = . TRUE . 

END IF 


C. . .End setup for wave simulator 

RETURN 

END 
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SUBROUTINE GETWAVE2 (METDEN, LATIN, APIN, ISEEDO, TOTDEN) 
IMPLICIT NONE 

INTEGER JSTEP, NRESET, ISEED, ISEEDO, APB IN, LATBIN, 

> OLDLATB IN 

REAL METDEN, LATIN, APIN, TOTDEN, DT, TIMESTEP (720 ) , 

> ELAPTIME, TMAX, OLDAP, WAVE 

LOGICAL RESET, PICKBINS, WARNINGS 

COMMON / WAVE / DT, TMAX, ISEED, RESET, PICKBINS, 

> APB IN, WARNINGS 

IF ( ELAPTIME. GT. TMAX .OR. (.NOT. RESET .AND . .NOT . PICKBINS 

> .AND. OLDAP .NE .APIN) ) RESET = .TRUE. 

IF (ABS (LATIN) .GT. 40 . 0) THEN 

LATBIN * 2 
ELSE 

LATBIN = 1 
END IF 

IF (LATB IN. NE. OLDLATB IN) RESET = .TRUE. 

IF (RESET) THEN 
NRESET = NRESET+1 

ISEED = ISEEDO +2 *NRESET ! this keeps seed odd 

OLDAP = APIN 
OLDLATB IN = LATBIN 
ELAPTIME = 0.0 
END IF 

C. . .Note that WAVES is called every time 

CALL WAVES 2 (APIN, LATIN, DT, TMAX, ISEED, RESET, 

> PICKBINS, APB IN, WARNINGS, WAVE) 

ELAPTIME = ELAPTIME +DT 

TOTDEN = (1.0 +WAVE ) * METDEN 

RETURN 

END 
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SUBROUTINE WAVES 2 (APIN, LATIN, DT, TMAX, ISEED , RESET, 

PICKBINS, APB IN, WARNINGS, WAVE) 


This subroutine /model is based on the analysis of AEC NATE data (see 
accompanying report by Michael P. Hickey, 1995) . 

It is valid only for inclinations greater than about 45 degrees 
and for approximately circular orbits between about 200 and 500 kra 
The time resolution of the original data and for this model is 
15 . 0 seconds . The perturbations output by this subroutine are used 
to modify the total mass density of the MET model and are intended 
for simulations related to atmospheric drag effects (e.g., guidance, 
navigation and control) . 


NPUTS : 


APIN 

LATIN 

DT 


TMAX 

ISEED 

RESET 

PICKBINS 

APB IN 

WARNINGS 


is the same ap value that is input to MET model; 
is latitude; 

is the next time step, which must be at least 
15 seconds (the time resolution of original data) . 
If the input value of DT is less than 15 s it is 
automatically reset to 15.0; 

is the maximum time duration (sec) of simulation. 

It is advised that TMAX not exceed 3 hours (the 
time resolution of the ap index) ; 

is a large odd integer for random number generator; 
tells subroutine to go back to setup mode if it is 
true. It should be reset when the ap bins varies ; 
tells subroutine to override APIN inputs and 
explicitly select required ap bin; 

is* the explicitly selected ap bin when PICKBINS is 
true ; 

tells subroutine to type all warnings that occur 
when either invalid input parameters are reset 
internally or when WAVES 2 should be recalled in 
setup mode. If WARNINGS is false, no warnings are 
typed to the screen, although the user should be 
aware that internal resetting may be occuring if 
the subroutine is being misused. 


OUTPUT: WAVE is the wave percentage amplitude expressed as a 

decimal . WAVE is a single dependant variable . 
Subroutine WAVES 2 is called every time step . 


C. . .Additional Subroutines 


Employed: GAUSSD and RANI (random 


# 


generators ) 
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IMPLICIT NONE 

INTEGER I SEED, LATBIN, APB IN, N, K, I, I GRID, JGRID 

REAL TMAX, 

> APIN, LATIN, DT, A, Al, A 2, WAVE, NOISE(721), PHI, 

> STNDEV, RATIO, ELTIME, PHI1, PHI2, 

> SCALE (2,2) , AR2 (8) 

LOGICAL FIXDSTEP, SETUP, RESET, PICKBINS, WARNINGS, 

> STOPFLAG 

2 . . . AE-C parameters: 

DATA AR2 / 1.1792, -0.24832, 1.2764, -0.33955, 

> 1.4367, -0.51164, 1.3461, -0.42737 / 

DATA SCALE / 1.385E-2, 2.093E-2, 1.171E-2, 2.440E-2 / 

DATA SETUP / . TRUE . / 

IF (RESET) SETUP = .TRUE. 

-• • -Check that DT and TMAX are both nonzero. Type fatal error messaae 

1.. . and terminate program execution if either one of them is zero. 

STOPFLAG = . FALSE . 

IF (DT.EQ.O.) THEN 
PRINT 1 

STOPFLAG = . TRUE . 

ELSE IF (TMAX.EQ.O.) THEN 
PRINT 2 

STOPFLAG = .TRUE. 

END IF 

IF (STOPFLAG) STOP 

— • • .Further check that DT is no smaller than 15 sec and that TMAX is no 
C. . .greater than 3 hrs . If they are, reset them to their nominal values 
IF (DT . LT .15.0) THEN 
DT = 15.0 

IF (WARNINGS) PRINT 3 
END IF 

IF (TMAX.GT. 1. 08E4) THEN 
TMAX = 1.08E4 
IF (WARNINGS) PRINT 4 
END IF 

C. . . Commence setting latitide & ap bins, white noise array on a constant 
C... 15-sec time grid. 

IF (SETUP) THEN 

SETUP = . FALSE . 

RESET = . FALSE . 

N = 1+TMAX/1S . 0 

C. . .Generate white noise of zero mean and unit standard deviation 
CALL GAUSSD (NOISE, N, I SEED) 

C. . .PICKBINS = true will override the input value of aD. 

IF ( .NOT. PICKBINS) THEN 
IF (APIN.LE. 7.0) THEN 
APB IN = 1 
ELSE 

APB IN = 2 
END IF 
END IF 

IF (ABS (LATIN) . LE . 40 . 0) THEN 

CIO 



o n 


c. 


LATBIN = 1 

PHI1 = AR2 (2*APBIN-1) 

PHI 2 = AR2 (2*APBIN) 

ELSE 

LATBIN = 2 

PHI1 = AR2 (4+2*APBIN-l) 

PHI 2 = AR2 (4+2*APBIN) 

END IF 

STNDEV = SCALE (LATBIN, APB IN) 
END IF 


.Wave model is valid on a fixed 15 second grid, 

Hmp steC5 can be input. We must therefore keep track or tne uocai 

leased ?im“ and a? each (possibly variable) step determine the 

.corresponding position on the fixed 1 .,f c ar ^ d wh" 1“ 

Mnt . p -u at T GRID is the position on the fixed 15 -sec gr.a, wn.io 

'.JGMD^'is the projected position calculated from tne elapsed time. 

.For first time: 

IF (ELTIME.EQ. 0 • 0) THEN 
I GRID = 1 
JGRID = 1 
A = NOISE (1) 

.and for all subsequent times: 


ELSE 

JGRID = 1+INT (ELTIME/15 . 0) 
DO WHILE (IGRID.LT. JGRID) 

I GRID = IGRID+1 
A = PHI1*A1 + PHI2*A2 + 
END DO 
END IF 
A2 = A1 
A1 = A 


NOISE (IGRID) 


Now make variance fit data by increasing s 
unity to STNDEV . 

WAVE = STNDEV*A 


andard deviation 


from 


2 

3 

4 

5 


IF (ELTIME .GT . TMAX) THEN 
PRINT 5 
WAVE =0.0 
END IF 

ELTIME = ELTIME+DT 


FORMAT ( ‘ 
FORMAT ( ' 
FORMAT ( ' 
FORMAT ( 
FORMAT ( 


FATAL ERROR: DT is zero in WAVES2' ) 

FATAL ERROR: TMAX is zero seconds o 

WARNING: DT has been increased to 15 s econcs ) 
WARNING • TMAX has been decreased to 3 nours ) 

Elaosed time > TMAX. WAVE now set to 0 
Recall WAVES 2 in setup mode ) 


0 ' , /, 


RETURN 

END 
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SUBROUTINE GAUSSD (X, N, ix) 

routine to produce vector of gaussian distributed random numbers 
in X(n) . The series in X(n) will have approximate zero mean and 
standard deviation of 1.0. Polar Method algorithm taken from 
p. 104, Seminumerical Algorithms, 

vol 2 of The Art of Computer Programming by D. E. Knuth, 196 9. 

ix is random number seed, i.e., large odd integer. 

implicit none 
integer i, ix, iy, n 
REAL X (N) , RANI 
real s, vl, v2 


1 = 0 
iy = ix 

do while (i .It. n) 

VI = 2.0*ranl(iy) - 1.0 
V2 = 2.0*ranl(iy) - 1.0 

S = VI* VI + V2*V2 

if (s .It. 1.0) then 

S = SQRT (-2.0 * ALOG ( S ) / S) 

1 = 1 + 1 
X(I) = VI * s 
1 = 1 + 1 

if (i .le. n) then 
X(I) = V2 * S 
end if 

end if 

end do 

ix = iy 
RETURN 

END 
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FUNCTION RANI (IDUM) 

3 Returns a uniform random deviate between 0.0 and 1.0. Set IDUM to any 
3 negative value to initialize or reinitialize the sequence. 

3 Copied from Numerical Recipes, W.H. Press, B.P. Flannery, 

3 S . A. Teukolsky & W.T. Vetterling, Cambridge Univ. Press, 1986, p.196 
IMPLICIT NONE 

INTEGER IFF, 1X1, 1X2, 1X3, J, Ml, M2, M3, IA1 , IA2 , IA3 , 

> IC1, IC2 , IC3 , IDUM 

REAL RM1 , RM2 , R, RANI 

DIMENSION R (97) 

PARAMETER (Ml = 259200, IA1-7141, IC1=54773, RM1-1./M1) 

PARAMETER (M2=134456, IA2 = 8121, IC2=28411, RM2 = 1./M2) 

PARAMETER (M3=243000, IA3-4561, IC3=51349) 

DATA IFF /0/ ! as above, initialize on first call even if 

! IDUM is not negative 
IF (IDUM.LT.O .OR. IFF.EQ.O) THEN 
IFF-1 

1X1 =MOD { I Cl - IDUM , Ml ) ! Seed the first routine, 

1X1 =MOD ( IA1 * 1X1 + I Cl , Ml ) 

1X2 =MOD ( 1X1 , M2 ) ! and use it to seed the second 

IXl=MOD ( IA1* IX1+IC1 , Ml ) 

1X3 =MOD ( 1X1 , M3 ) ! and third routines. 

DO J = 1, 97 ! Fill the table with sequential 

IXl=MOD ( IA1*IX1+IC1 , Ml ) ! uniform deviates generated by 

IX2=MOD (IA2 *1X2+1 C2, M2) ! the first two routines 

R (J) = (FLOAT (1X1) +FLOAT (1X2) *RM2) *RM1 ! Low and high- 

END DO ! order pieces combined here . 

IDUM-1 
END IF 

IXl=MOD (IA1*IX1+IC1,M1) ! Except when initializing, this is 

IX2=MOD (IA2 *1X2+1 C2, M2) ! where we start. Generate the next 

IX3=MOD (IA3 *1X3 + 1 C3 , M3) ! number for each sequence. 

j=l+ (97*1X3) /M3 ! Use the third sequence to get an 

! integer between 1 and 97 
IF ( J . GT . 97 .OR. J.LT.l) PAUSE 

RAN1=R ( J) ! return that table entry, 

R(J) = (FLOAT (1X1) +FLOAT( 1X2) *RM2) *RM1 ! and refill it. 

RETURN 

END 
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APPENDIX D 


USER'S SOFTWARE IMPLMENTATION GUIDE 

AND 

PROGRAM LISTINGS 

THE MET MODEL 
STATISTICAL ANALYSIS MODE 
(MET-SAM) 



USER’S SOFTWARE IMPLEMENTATION GUIDE 


To use the statistical analysis mode(SAM) of the MET model select the appro- 
priate exospheric temperature from the following tables and the altitude for which you 
desire the density to be computed. Enter these two values into the subroutine program 
code. To obtain an exospheric temperature value at a specific percentile level that is not in 
the table, linearly interpolate between the two nearest values. PFD is the percent 
frequency distribution. 


Table 1. Cumulative Percent Frequency Distribution of Minimum Exospheric 
Temperature When 13-Month Smoothed 10.7-cm Solar Flux Is Between 66 and 246 


PFD 

CPF 


PFD 

CPF 


PFD 

CPF 


0.0256 

0.02 

567.00 

2.59 

75.21 

987.00 

0.16 

9937 

1407.00 

0.29 

0.31 

587.00 

2.54 

77.74 

1007.00 

0.17 

9934 

1427.00 

1.34 

1.66 

607.00 

2.29 

80.04 

1027.00 

0.13 

99.67 

1447.00 

2.70 

4.36 

627.00 

2.30 

82.34 

1047.00 

0.09 

99.76 

1467.00 

4.33 

8.69 

647.00 

2.24 

84.57 

1067.00 

0.06 

99.82 

1487.00 

5.87 

14.56 

667.00 

1.97 

86.54 

1087.00 

0.05 

99.88 

1507.00 

6.11 

20.68 

687.00 

1.80 

88.34 

1107.00 

0.05 

99.92 

1527.00 

5.40 

26.08 

708.00 

1.60 

89.93 

1127.00 

0.03 

99.95 

1547.00 

5.01 

31.09 

727.00 

1.49 

91.43 

1147.00 

0.02 

99.97 

1567.00 

4.79 

35.88 

747.00 

1.34 

92.77 

1167.00 

0.01 

99.98 

1587.00 

4.43 

40.31 

767.00 

1.18 

93.95 

1187.00 

0.01 

99.98 

1607.00 

4.05 

44.36 

787.00 

1.11 

95.06 

1207.00 

0.01 

99.99 

1627.00 

3.61 

47.98 

807.00 

0.94 

96.00 

1227.00 

0.00 

99.99 

1647.00 

3.34 

51.31 

827.00 

0.73 

96.73 

1247.00 

0.01 

100.00 

1667.00 

3.27 

54.59 

847.00 

0.60 

97.33 

1267.00 

0.00 

100.00 

1687.00 

3.38 

57.96 

867.00 

0.50 

97.82 

1287.00 

0.00 

100.00 

1707.00 

3.25 

61.21 

887.00 

0.41 

98.23 

1307.00 

0.00 

100.00 

1727.00 

3.09 

64.30 

907.00 

0.36 

98.59 

1327.00 

0.00 

100.00 

1747.00 

3.05 

67.35 

927.00 

0.25 

98.83 

1347.00 

0.00 

100.00 

1767.00 

2.61 

69.96 

947.00 

0.18 

99.02 

1367.00 


2.65 

72.61 

967.00 

0.19 

99.21 

1387.00 
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Table 2. Cumulative Percent Frequency Distribution of Mean Exospheric 
Temperature When 13-Month Smoothed 10.7-cm Solar Flux Is Between 66 and 246 
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Table 3. Cumulative Percent Frequency Distribution of Maximum Exospheric 
Temperature When 13-Month Smoothed 10.7-cm Solar Flux Is Between 66 and 246 
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Table 4. Cumulative Percent Frequency Distribution of Minimum Exospheric 
Temperature When 13-Month Smoothed 10.7-cm Solar Flux Is Between 66 and 102 


PFD 

CPF 


PFD 

CPF 


PFD 

CPF 


"004 

0.04 

567.00 

4.23 

92.67 

787.00 

0.02 

99.98 

1007.00 

0.74 

0.78 

587.00 

2.85 

95.53 

807.00 

0.01 

99.98 

1027.00 

3.39 

4.17 

607.00 

1.72 

97.24 

827.00 

0.00 

99.98 

1047.00 

6.80 

10.98 

627.00 

1.13 

98.37 

847.00 

0.01 

99.99 

1067.00 

10.84 

21.82 

647.00 

0.72 

99.09 

867.00 

0.01 

100.00 

1087.00 

14.65 

36.47 

667.00 

0.35 

99.44 

887.00 

0.00 

100.00 

1107.00 

14.84 

51.31 

687.00 

0.25 

99.69 

907.00 

0.00 

100.00 

1127.00 

12.25 

63.56 

707.00 

0.13 

99.82 

927.00 

0.00 

100.00 

1147.00 

10.17 

73.73 

727.00 

0.07 

99.88 

947.00 

0.00 

100.00 

1167.00 

8.36 

82.09 

747.00 

0.06 

99.94 

967.00 

0.00 

100.00 

1173.00 | 

6.36 

88.44 

767.00 

0.01 

99.96 

987.00 

I 


Table 5. Cumulative Percent Frequency Distribution of Mean Exospheric 
Temperature When 13-Month Smoothed 10.7-cm Solar Flux Is Between 66 and 102 


| PFD 

CDF 


PFD 

CDF 


PFD 

CDF 

1 

0.04 

0.04 

660.00 

5.03 

90.32 

880.00 

0.03 

99.96 

1100.00 

0.62 

0.66 

680.00 

3.48 

93.80 

900.00 

0.01 

99.97 

1120.00 

2.81 

3.47 

700.00 

2.25 

96.05 

920.00 

0.01 

99.98 

1140.00 

5.73 

9.20 

720.00 

1.42 

97.47 

940.00 

0.01 

99.99 

1160.00 

9.33 

18.54 

740.00 

1.04 

98.50 

960.00 

0.00 

99.99 

1180.00 

13.23 

31.76 

760.00 

0.66 

99.16 

980.00 

0.01 

100.00 

1200.00 

14.35 

46.12 

780.00 ' 

0.29 

99.45 

1000.00 

0.00 

100.00 

1220.00 

12.76 

58.87 

800.00 

0.22 

99.67 

1020.00 

0.00 

100.00 

1240.00 

10.80 

69.67 

820.00 

0.14 

99.82 

1040.00 

0.00 

100.00 

1260.00 

8.82 

78.49 

840.00 

0.08 

99.90 

1060.00 

0.00 

100.00 

1278.00 

6.80 

85.29 

860.00 

0.04 

99.94 

1080.00 

1 
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Table 6. Cumulative Percent Frequency Distribution of Maximum Exospheric 
Temperature When 13-Month Smoothed 10.7-cm Solar Flux is Between 66 and 102 



Table 7. Cumulative Percent Frequency Distribution of Minimum Exospheric 
Temperature When 13-Month Smoothed 10.7 Solar Flux is Between 102 and 138 
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Table 8. Cumulative Percent Frequency Distribution of Mean Exospheric 
Temperature When 13-Month Smoothed 10.7-cm Solar Flux Is Between 102 and 138 


Iffd" 

CPF 


PFD 

CPF 


PFD 

CPF 


1 0.03 

0.03 

718.00 

4.99 

85.55 

1038.00 

0.00 

99.93 

1358.00 

0.08 

0.11 

738.00 

3.67 

89.22 

1058.00 

0.01 

99.94 

1378.00 

0.22 

0.32 

758.00 

2.96 

92.18 

1078.00 

0.01 

99.95 

1398.00 

0.36 

0.68 

778.00 

2.42 

94.60 

1098.00 

0.01 

99.96 

1418.00 

1.18 

1.86 

798.00 

1.76 

96.37 

1118.00 

0.00 

99.96 

1438.00 

2.53 

4.39 

818.00 

1.22 

97.59 

1138.00 

0.01 

99.98 

1458.00 

4.34 

8.74 

838.00 

0.80 

98.39 

1158.00 

0.00 

99.98 

1478.00 

6.21 

14.95 

858.00 

0.48 

98.87 

1178.00 

0.01 

99.99 

1498.00 

7.51 

22.46 

878.00 

0.36 

99.23 

1198.00 

0.00 

99.99 

1518.00 

8.81 

31.27 

898.00 

0.28 

99.50 

1218.00 

0.00 

99.99 

1538.00 

9.14 

40.41 

918.00 

0.13 

99.64 

1238.00 

0.00 

99.99 

1558.00 

935 

49.76 

938.00 

0.08 

99.72 

1258.00 

0.00 

99.99 

1578.00 

8.55 

58.31 

958.00 

0.10 

99.82 

1278.00 

0.01 

100.00 

1598.00 

8.08 

66.39 

978.00 

0.06 

99.88 

1298.00 

0.00 

100.00 

1608.00 

7.69 

74.08 

998.00 

0.04 

99.92 

1318.00 


6.48 

80.56 

1018.00 

0.00 

99.93 

1338.00 
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Table 9. Cumulative Percent Frequency Distribution of Maximum Exospheric 
Temperature When 13-Month Smoothed 10.7-cm Solar Flux Is Between 102 and 138 


PFD 

CPF 


PFD 

CPF 


PFD 

CPF 


0.03 

0.03 

820.00 

6.28 

79.22 

1140.00 

0.04 

99.91 

1460.00 

0.08 

0.11 

840.00 

4.97 

84.19 

1160.00 

0.02 

99.93 

1480.00 

| 0.18 

0.28 

860.00 

3.86 

88.05 

1180.00 

0.00 

99.93 

1500.00 

0.24 

0.52 

880.00 

2.65 

90.70 

1200.00 

0.01 

99.94 

1520.00 

0.64 

1.16 

900.00 

2.44 

93.13 

1220.00 

0.01 

99.95 

1540.00 

2.00 

3.16 

920.00 

2.13 

95.26 

1240.00 

0.00 

99.96 

1560.00 

3.17 

6.33 

940.00 

1.57 

96.83 

1260.00 

0.01 

99.96 

1580.00 

4.87 

11.21 

960.00 

1.00 

97.83 

1280.00 

0.01 

99.98 

1600.00 

6.57 

17.78 

980.00 

0.65 

98.48 

1300.00 

0.00 

99.98 

1620.00 

7.55 

25.32 

1000.00 

0.40 

98.88 

1320.00 

0.01 

99.99 

1640.00 

8.26 

33.58 

1020.00 

0.33 

99.21 

1340.00 

0.00 

99.99 

1660.00 

8.66 

42.24 

1040.00 

0.24 

99.45 

1360.00 

0.00 

99.99 

1680.00 

8.85 

51.09 

1060.00 

0.17 

99.63 

1380.00 

0.00 

99.99 

1700.00 

7.96 

59.05 

1080.00 

0.08 

99.71 

1400.00 

0.00 

99.99 

1720.00 

7.07 

66.12 

1100.00 

0.09 

99.80 

1420.00 

0.01 

100.00 

1740.00 

6.82 

72.94 

1120.00 

0.07 

99.87 

1440.00 

0.00 

100.00 

1749.00 
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Table 10. Cumulative Percent Frequency Distribution of Minimum Exospheric 
Temperature When 13-Month Smoothed 10.7-cm Solar Flux Is Between 138 and 174 


PFD 

CPF 


PFD 

CPF 


PFD 

CPF 


0.08 

0.08 

714.00 

7.66 

70.07 

994.00 

0.12 

99.64 

1274.00 

031 

0.39 

734.00 

6.70 

76.77 

1014.00 

0.14 

99.79 

1294.00 

1.14 

1.53 

754.00 

6.06 

82.84 

1034.00 

0.05 

99.84 

1314.00 

2.07 

3.60 

774.00 

4.77 

87.60 

1054.00 

0.07 

99.91 

1334.00 

2.24 

5.84 

794.00 

3.79 

91.40 

1074.00 

0.03 

99.94 

1354.00 

3.02 

8.87 

814.00 

2.47 

93.87 

1094.00 

0.02 

99.96 

1374.00 

4.16 

13.03 

834.00 

1.80 

95.66 

1114.00 

0.01 

99.97 

1394.00 

5.32 

18.35 

854.00 

1.39 

97.06 

1134.00 

0.01 

99.98 

1414.00 

5.94 

24.29 

874.00 

0.83 

97.88 

1154.00 

0.00 

99.99 

1434.00 

6.58 

30.87 

894.00 

0.55 

98.43 

1174.00 

0.01 

100.00 

1454.00 

8.15 

39.02 

914.00 

0.40 

98.82 

1194.00 

0.00 

100.00 

1474.00 

836 

47.28 

934.00 

0.29 

99.12 

1214.00 

0.00 

100.00 

1494.00 

7.69 

54.97 

954.00 

0.19 

99.31 

1234.00 


1 7.44 

62.41 

974.00 

031 

99.52 

1254.00 
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Table 11. Cumulative Percent Frequency Distribution of Mean Exospheric 
Temperature When 13-Month Smoothed 10.7-cm Solar Flux Is Between 138 and 174 


PFD 

CPF 


PFD 

CPF 


PFD 

CPF 


0.08 

0.08 

836.00 

7.00 

69.56 

1136.00 

0.13 

99.59 

1436.00 

031 

039 

856.00 

6.51 

76.07 

1156.00 

0.14 

99.73 

1456.00 

0.95 

1.34 

876.00 

5.63 

81.71 

1176.00 

0.09 

99.82 

1476.00 

1.67 

3.01 

896.00 

4.70 

86.41 

1196.00 

0.04 

99.86 

1496.00 

2.20 

5.21 

916.00 

3.89 

90.30 

1216.00 

0.06 

99.92 

1516.00 

2.54 

7.74 

936.00 

2.59 

92.90 

1236.00 

0.03 

99.96 

1536.00 

335 

11.09 

956.00 

1.91 

94.80 

1256.00 

0.02 

99.98 

1556.00 

4.50 

15.60 

976.00 

1.46 

96.27 

1276.00 

0.00 

99.99 

1576.00 

5.38 

20.98 

996.00 

1.02 

97.29 

1296.00 

0.00 

99.99 

1596.00 

5.61 

26.58 

1016.00 

0.70 

97.99 

1316.00 

0.01 

100.00 

1616.00 

6.61 

33.19 

1036.00 

0.52 

98.51 

1336.00 

0.00 

100.00 

1636.00 

7.93 

41.12 

1056.00 

034 

98.84 

1356.00 

0.00 

100.00 

1656.00 

736 

48.68 

1076.00 

0.28 

99.12 

1376.00 

0.00 

100.00 

1659.00 

7.03 

55.71 

1096.00 

0.17 

99.29 

1396.00 


6.85 

62.56 

1116.00 

0.17 

99.46 

1416.00 
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Table 12. Cumulative Percent Frequency Distribution of Maximum Exospheric 
Temperature When 13-Month Smoothed 10.7-cm Solar Flux Is Between 138 and 174 


PFD 

CPF 


PFD 

CPF 


PFD 

CPF 

1 

0.08 

0.08 

958.00 

6.35 

62.90 

1258.00 

0.21 

9930 

1558.00 1 

0.33 

0.41 

978.00 

6.22 

69.11 

1278.00 

0.14 

99.43 

1578.00 1 

0.83 

1.24 

998.00 

6.11 

75.22 

1298.00 

0.15 

99.58 

1598.00 

1.40 

2.64 

1018.00 

5.47 

80.69 

1318.00 

0.07 

99.65 

1618.00 

2.04 

4.68 

1038.00 

4.48 

85.17 

1338.00 

0.11 

99.76 

1638.00 

2.13 

6.81 

1058.00 

3.76 

88.93 

1358.00 

0.08 

99.84 

1658.00 

2.98 

9.79 

1078.00 

2.97 

91.90 

1378.00 

0.05 

99.90 

1678.00 

3.73 

13.51 

1098.00 

1.95 

93.85 

1398.00 

0.04 

99.94 

1698.00 

4.54 

18.05 

1118.00 

1.54 

95.39 

1418.00 

0.03 

99.97 

1718.00 

5.28 

23.33 

1138.00 

1.21 

96.60 

1438.00 

0.02 

99.99 

1738.00 

5.43 

28.75 

1158.00 

0.85 

97.45 

1458.00 

0.00 

99.99 

1758.00 

6.70 

35.46 

1178.00 

0.61 

98.06 

1478.00 

0.01 

■aa 

1778.00 

7.38 

42.84 

1198.00 

0.45 

9852 

1498.00 

0.00 

100.00 

1798.00 

6.99 

49.83 

1218.00 

0.33 

98.84 

1518.00 

0.00 

100.00 

1818.00 

6.71 

56.34 

1238.00 

0.25 

99.09 

1538.00 

0.00 

100.00 

1825.00 
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Table 13. Cumulative Percent Frequency Distribution of Minimum Exospheric 
Temperature When 13-Month Smoothed 10.7-cm Solar Flux Is Between 174 and 210 


PFD 

CPF 


PFD 

CPF 


PFD 

CPF 


0.02 

0.02 

764.00 

6.99 

55.04 

1084.00 

0.26 

99.18 

1404.00 

0.11 

0.13 

784.00 

6.60 

61.64 

1104.00 

0.18 

9936 

1424.00 

0.12 

0.25 

804.00 

5.87 

67.51 

1124.00 

0.24 

9939 

1444.00 

0.29 

0.54 

824.00 

5.81 

73.32 

1144.00 

0.10 

99.70 

1464.00 

0.72 

1.27 

844.00 

4.99 

78.32 

1164.00 

0.08 

99.78 

1484.00 

1.53 

2.79 

864.00 

4.16 

82.48 

1184.00 

0.07 

99.85 

1504.00 

2.79 

5.58 

884.00 

3.99 

86.47 

1204.00 

0.04 

99.89 

1524.00 

3.25 

8.83 

904.00 

3.48 

89.94 

1224.00 

0.06 

99.96 

1544.00 

3.45 

12.28 

924.00 

2.65 

92.90 

1244.00 

0.03 

99.98 

1564.00 

3.12 

15.39 

944.00 

1.94 

94.54 

1264.00 

0.01 

99.99 

1584.00 

4.23 

19.62 

964.00 

1.45 

95.99 

1284.00 

0.00 

99.99 

1604.00 

4.78 

24.40 

984.00 

1.19 

97.18 

1304.00 

0.00 

99.99 

1624.00 

534 

29.74 

1004.00 

0.70 

97.89 

1324.00 

0.00 

99.99 

1644.00 

5.65 

35.40 

1024.00 

0.52 

98.41 

1344.00 

0.01 

99.99 

1664.00 

6.17 

41.56 

1044.00 

0.27 

98.67 

1364.00 

0.01 

100.00 

1684.00 

6.49 

48.05 

1064.00 

0.24 

98.92 

1384.00 

0.00 

100.00 

1696.00 
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Table 14. Cumulative Percent Frequency Distribution of Mean Exospheric 
Temperature When 13-Month Smoothed 10.7-cm Solar Flux Is Between 174 and 210 


PFD 

CPF 


PFD 

CPF 


PFD 

CPF 


0.03 

0.03 

894.00 

6.29 

52.11 

1234.00 

0.15 

98.93 

1574.00 

0.10 

0.13 

914.00 

6.23 

58.35 

1254.00 

0.22 

99.16 

1594.00 

0.10 

0.23 

934.00 

5.75 

64.10 

1274.00 

0.22 

99.38 

1614.00 

0.19 

0.42 

954.00 

5.40 

69.50 

1294.00 

0.21 

99.58 

1634.00 

0.41 

0.83 

974.00 

4.91 

74.41 

1314.00 

0.11 

99.70 

1654.00 

1.08 

1.91 

994.00 

438 

78.79 

1334.00 

0.08 

99.77 

1674.00 

1.99 

3.90 

1014.00 

3.98 

82.76 

1354.00 

0.04 

99.82 

1694.00 

0.03 

0.03 

894.00 

6.29 

52.11 

1234.00 

0.15 

98.93 

1574.00 

0.10 

0.13 

914.00 

6.23 

5835 

1254.00 

0.22 

99.16 

1594.00 

0.10 

0.23 

934.00 

5.75 

64.10 

1274.00 

0.22 

99.38 

1614.00 

0.19 

0.42 

954.00 

5.40 

69.50 

1294.00 

0.21 

9938 

1634.00 

0.41 

0.83 

974.00 

4.91 

74.41 

1314.00 

0.11 

99.70 

1654.00 

1.08 

1.91 

994.00 

4.38 

78.79 

1334.00 

0.08 

99.77 

1674.00 

1.99 

3.90 

1014.00 

3.98 

82.76 

1354.00 

0.04 

99.82 

1694.00 

2.98 

6.88 

1034.00 

3.48 

86.24 

1374.00 

0.06 

99.88 

1714.00 

3.20 

10.08 

1054.00 

3.33 

89.57 

1394.00 

0.07 

99.95 

1734.00 

3.06 

13.14 

1074.00 

2.67 

92.24 

1414.00 

0.02 

99.97 

1754.00 

2.87 

16.00 

1094.00 ! 

1.77 

94.02 

1434.00 

0.02 

99.99 

1774.00 

3.53 

19.53 

1114.00 

1.43 

95.45 

1454.00 

0.00 

99.99 

1794.00 

4.59 

24.12 

1134.00 

1.19 

96.64 

1474.00 

0.00 

99.99 

1814.00 

5.06 

29.19 

1154.00 

0.90 

97.54 

1494.00 

0.01 

99.99 

1834.00 

5.11 

34.30 

1174.00 

0.64 

98.18 

1514.00 

0.00 

99.99 

1854.00 

5.84 

40.14 

1194.00 

038 

98.56 

1534.00 

0.01 

100.00 

1874.00 

5.69 

45.83 

1214.00 

0.23 

98.79 

1554.00 

0.00 

— 

1876.00 
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Table 15. Cumulative Percent Frequency Distribution of Maximum Exospheric 
Temperature When 13-Month Smoothed 10.7-cm Solar Flux Is Between 174 and 210 


PFD 

CPF 


PFD 

CPF 


PFD 

CPF 

| 

0.03 

0.03 

1024.00 

5.63 

49.93 

1384.00 

0.16 

98.82 

1744.00 | 

0.10 

0.12 

1044.00 

5.63 

55.56 

1404.00 

0.13 

98.95 

1764.00 1 

0.10 

0.23 

1064.00 

5.45 

61.01 

1424.00 

0.20 

99.15 

1784.00 1 

0.14 

0.36 

1084.00 

5.34 

66.35 

1444.00 

0.21 

99.36 

1804.00 

0.25 

0.61 

1104.00 

4.76 

71.11 

1464.00 

0.19 

99.55 

1824.00 

0.71 

1.33 

1124.00 

4.26 

75.37 

1484.00 

0.15 

99.70 

1844.00 

1.44 

2.77 

1144.00 

3.81 

79.18 

1504.00 

0.06 

99.76 

1864.00 

2.43 

5.20 

1164.00 

3.67 

82.85 

1524.00 

0.05 

99.81 

1884.00 

2.91 

8.11 

1184.00 

3.27 

86.13 

1544.00 

0.05 

99.85 

1904.00 

3.05 

11.16 

1204.00 

3.03 

89.16 

1564.00 

0.06 

99.91 

1924.00 

2.73 

13.39 

1224.00 

2.55 

91.71 

1584.00 

0.04 

99.95 

1944.00 

2.73 

16.62 

1244.00 

1.93 

93.64 

1604.00 

0.03 

99.98 

1964.00 

3.09 

19.71 

1264.00 

1.34 

94.98 

1624.00 

0.01 

99.99 

1984.00 

4.34 

24.05 

1284.00 

1.18 

96.16 

1644.00 

0.01 

99.99 

2004.00 

4.66 

28.72 

1304.00 

0.95 

97.11 

1664.00 

0.00 

99.99 

2024.00 

5.03 

33.75 

1324.00 

0.78 

97.89 

1684.00 

0.01 

100.00 

2044.00 

5.15 

38.90 

1344.00 

0.49 

9837 

1704.00 

0.00 

100.00 

2055.00 

5.40 

44.30 

1364.00 

0.28 

98.65 

1724.00 

1 
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Table 16. Cumulative Percent Frequency Distribution of Minimum Exospheric 
Temperature When 13-Month Smoothed 10.7-cm Solar Flux Is Between 210 and 246 


Ipfd 

CPF 


PFD 

CPF 


PFD 

CPF 

1 

| 0.03 

0.03 

924.00 

5.50 

61.90 

1224.00 

0.50 

98.89 

1524.00 

| 034 

0.37 

944.00 

5.07 

66.97 

1244.00 

034 

99.23 

1544.00 

I 0.69 

1.06 

964.00 

4.41 

71.39 

1264.00 

037 

99.50 

1564.00 

1.26 

2.31 

984.00 

3.89 

75.27 

1284.00 

0.16 

99.66 

1584.00 J 

1.30 

3.61 

1004.00 

3.44 

78.71 

1304.00 

0.11 

99.77 

1604.00 

1.94 

5.56 

1024.00 

3.67 

82.89 

1324.00 

0.04 

99.81 

1624.00 

3.43 

8.99 

1044.00 

2.80 

85.19 

1344.00 

0.07 

99.89 

1644.00 

4.93 

13.91 

1064.00 

2.20 

87.39 

1364.00 

0.09 

99.97 

1664.00 

5.26 

19.17 

1084.00 

2.37 

89.76 

1384.00 

0.01 

99.99 

1684.00 

5.44 

24.61 

1104.00 

1.91 

91.67 

1404.00 

0.00 

99.99 

1704.00 

5.50 

30.11 

1024.00 

2.17 

93.84 

1424.00 

0.00 

99.99 

1724.00 

I 6 * 21 

36.33 

1144.00 

1.59 

95.43 

1444.00 

0.00 

99.99 

1744.00 

1 6.84 

43.17 

1164.00 

1.27 

96.70 

1464.00 

0.01 

100.00 

1764.00 

1 6.83 

50.00 

1184.00 

0.86 

97.56 

1484.00 

0.00 

100.00 

1767.00 

[6.40 

56.40 

1204.00 

0.83 

98.39 

1504.00 

J 
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Table 17. Cumulative Percent Frequency Distribution of Mean Exospheric 
Temperature When 13-Month Smoothed 10.7-cm Solar Flux Is Between 210 and 246 


1 PFD 

CPF 


PFD 

CPF 


PFD 

CPF 


0.03 

0.03 

1084.00 

4.06 

62.06 

1404.00 

0.47 

98.74 

1724.00 

030 

0.33 

1104.00 

4.49 

66.54 

1424.00 

0.41 

99.16 

1744.00 

037 

0.90 

1124.00 

4.24 

70.79 

1444.00 

031 

99.37 

1764.00 

1.21 

2.11 

1144.00 

3.69 

74.47 

1464.00 

0.20 

99.57 

1784.00 

1.06 

3.17 

1164.00 

3.34 

77.81 

1484.00 

0.17 

99.74 

1804.00 

1.60 

4.77 

1184.00 

3.20 

81.01 

1504.00 

0.06 

99.80 

1824.00 

2.50 

7.27 

1204.00 

2.86 

83.87 

1524.00 

0.06 

99.86 

1844.00 

4.53 

11.80 

1224.00 

2.03 

85.90 

1544.00 

0.10 

99.96 

1864.00 

5.49 

17.29 

1244.00 

2.01 

87.91 

1564.00 

0.01 

99.97 

1884.00 

4.74 

22.03 

1264.00 

2.24 

90.16 

1584.00 

0.01 

99.99 

1904.00 

5.56 

27.59 

1284.00 

1.97 

92.13 

1604.00 

0.00 

99.99 

1924.00 

5.40 

32.99 

1304.00 

1.81 

93.94 

1624.00 

0.00 

99.99 

1944.00 

6.11 

39.10 

1324.00 

1.64 

95.59 

1644.00 

0.01 

■n 

1964.00 

6.56 

45.66 

1344.00 

0.96 

96.54 

1664.00 

0.00 

100.00 

1975.00 

634 

52.00 

1364.00 

0.90 

97.44 

1684.00 


6.00 

58.00 

1384.00 

0.83 

98.27 

1704.00 
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Table 18. Cumulative Percent Frequency Distribution of Maximum Exospheric 
Temperature When 13-Month Smoothed 10.7-cm Solar Flux Is Between 210 and 246 


1 PFD 

CPF 


PFD 

CPF 


PFD 

CPF 

| 

I 0.03 

0.03 

1243.00 

3.89 

62.19 

1583.00 

0.64 

98.57 

1923.00 | 

0.21 

0.24 

1263.00 

3.64 

65.83 

1603.00 

0.40 

98.97 

1943.00 

0.56 

0.80 

1283.00 

4.06 

69.89 

1623.00 

036 

99.33 

1963.00 

0.97 

1.77 

1303.00 

3.60 

73.49 

1643.00 

0.09 

99.41 

1983.00 

1.11 

2.89 

1323.00 

3.31 

76.80 

1663.00 

0.19 

99.60 

2003.00 

1.39 

4.27 

1343.00 

3.01 

79.81 

1683.00 

0.14 

99.74 

2023.00 

1.94 

6.21 

1363.00 

2.49 

82.30 

1703.00 

0.13 

99.87 

2043.00 

3.39 

9.60 

1383.00 

2.36 

84.66 

1723.00 

0.06 

99.93 

2063.00 

5.44 

15.04 

1403.00 

1.97 

86.63 

1743.00 

0.03 

99.96 

2083.00 

4.86 

19.90 

1423.00 

1.71 

88.34 

1763.00 

0.03 

99.99 

2103.00 

5.19 

25.09 

1443.00 

2.19 

90.53 

1783.00 

0.00 

99.99 

2123.00 

4.67 

29.76 

1463.00 

1.89 

92.41 

1803.00 

0.00 

99.99 

2143.00 

1 5.91 

35.67 

1483.00 

1.41 

93.83 

1823.00 

0.00 

99.99 

2163.00 

— 

5.69 

41.36 

1503.00 

1.70 

95.53 

1843.00 

0.01 

100.00 

2183.00 

6.14 

47.50 

1523.00 

0.94 

96.47 

1863.00 

1 

5.89 

53.39 

1543.00 

0.80 

97.27 

1883.00 

I 

4.91 

58.30 

1563.00 

0.66 

97.93 

1903.00 
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PROGRAM LISTINGS 


SUBROUTINE J70MET (TE, Z, TZ, DENS, DL) 

C..J70MET developed from MET2 by Mike Hickey, Physitron Consultant 
C. INPUTS: 

C.. TE ~ exospheric temperature 

C.. Z — altitude 

C.. 

C.. OUTPUTS: 

C.. TZ — temperature at altitude Z 
C.. DENS -- total density at altitude Z 
C.. DL — loglO total density at altitude Z 

PARAMETER RGAS = 8.31432E3 ! J/kmol-K 

PARAMETER BFH = 440.0 

C Calcultions performed for only one latitude, one longitude 
C and one altitude 

CALL JAC (Z, TE, TZ, EM, DENS, DL) 

RETURN 

END 
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SUBROUTINE JAC (Z, T, TZ, EM, DENS, DL) 

C Subroutine 'JAC' calculates the temperature TZ, the total density DENS 
C and its logarithm DL, the mean molecular weight EM the individual 
C specie number densities for N, 02, 0, A, HE and H (each preceded with 
C an 'A') at altitude Z given the exospheric temperature T. 

C This subroutine uses the subroutine 'GAUSS' and the function 
C subprograms 'TEMP' and MOL_WT'. 

C 

C Originally Rewritten by Mike Hickey, USRA 

DIMENSION ALPHA(6), EI(6), DI(6), DIT(6) 

REAL *4 MOL WT 

PARAMETER AV = 6.02257E23 
PARAMETER QN = 78110 
PARAMETER Q02 =.20955 
PARAMETER QA =.009343 
PARAMETER QHE = 1.289E-5 
PARAMETER RGAS = 8.3 1432 
PARAMETER PI =3.14159265 
PARAMETER TO =183. 

GRAVITY (ALTITUDE) = 9.80665 / (( 1. + ALTITUDE / 6.356766E3)**2) 

DATA ALPHA / 0.0, 0.0, 0.0, 0.0, -.380, 0.0 / 

DATA El/ 28.0134, 31.9988, 15.9994, 39.948, 4.0026, 1.00797/ 

TX = 444.3807 + .02385 * T - 392.8292 * EXP (-.0021357 * T) 

A2 = 2. * (T-TX) / PI 
TX_T0 = TX - TO 
Tl= 1.9 * TX_T0 / 35. 

T3 = -1.7 * TX_T0 / (35.**3) 

T4 = -0.8 * TX_T0 / (3 5 * *4) 

TZ = TEMP (Z, TX, Tl, T3, T4, A2) 

C SECTION 1 


A = 90. 

D = AMIN1 (Z, 105.) 
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C Integrate gM/T from 90 to minimum of Z or 105 km 

CALL GAUSS (A, D, 1, R, TX, Tl, T3, T4, A2) 

C The number 2. 1926E-8 = density x temperature/mean molecular weight at 90 km. 
EM = MOL_WT (D) 

TD = TEMP (D, TX, Tl, T3, T4, A2) 

DENS = 2. 1926E-8*EM*EXP(-R/RGAS)/TD 
FACTOR = AV*DENS 
PAR = FACTOR /EM 
FACTOR = FACTOR / 28.96 

C For altitudes below and at 105 km calculate the individual specie number 
C densities from the mean molecular weight and total density. 

IF (Z. LE. 105) THEN 
DL = ALOG10 (DENS) 

AN = ALOG10 (QN*F ACTOR) 

AA = ALOG10 (QA*F ACTOR) 

AHE = ALOG10 (QHE*FACTOR) 

AO = ALOG10 (2.*PAR*(1.-EM/ 28.96)) 

A02 = ALOG10 (PAR*(EM*(1.+Q02) / 28.96-1.)) 

AH = 0. 

C Return to calling program 
RETURN 
END IF 

C SECTION 2 : This section is only performed for altitudes above 105 km 
C 

C Note that having reached this section means that D in section 1 is 105 km. 

C Calculate individual specie number densities from the total density and mean 
C molecular weight at 105 km altitude. 

DI(1) = QN*FACTOR 

DI(2) = PAR*(EM*(1 .+Q02) / 28.96-1.) 

DI(3) = 2.*PAR*(1.- EM / 28.96) 

DI(4) = QA*FACTOR 
DI(5) = QHE*F ACTOR 

C Integrate g/T from 105 km to Z km :- 

CALL GAUSS (D, Z, 2, R, TX, Tl, T3, T4, A2) 

DO 41 1= 1, 5 

DIT(I) = DI(I)*(TD/TZ)**(1.+ALPHA(I))*EXP(-EI(I)*R/RGAS) 
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IF (DIT(I). LE. 0.) DIT(I) = l.E-6 
41 CONTINUE 

C This section calculates atomic hydrogen densities above 500 km altitude. 
C Below this altitude, H densities are set to 10-6. 

C SECTION 3 


IF (Z GT. 500.) THEN 
A1 = 500. 

S = TEMP (Al, TX, Tl, T3, T4, A2) 

DI(6) = 10. **(73. 13 - 39.4*ALOG10 (S) + 5.5*ALOG10(S)*ALOG10(S)) 
CALL GAUSS (Al, Z, 7, R, TX, Tl, T3, T4, A2) 

DIT(6) = DI(6) * ( S/TZ) *EXP (-EI(6)*R / RGAS) 

ELSE 

DIT (6)= 1.0 
END IF 

C For altitudes greater than 105 km, calculate total density and mean 
C molecular weight from individual specie number densities. 

DENSO 
DO 42 1= 1, 6 

DENS = DENS + EI(I)*DIT(I) / AV 
42 CONTINUE 

EM = DENS*AV / (DIT(1)+DIT(2)+DIT(3)+DIT(4)+DIT(5)+DIT(6)) 

DL = ALOG10 (DENS) 

AN = ALOG10(DIT(1)) 

A02 = ALOG10(DIT(2)) 

AO = ALOG10(DIT(3)) 

AA = ALOG10(DIT(4)) 

AHE = ALOG10(DIT(5)) 

AH = ALOG10(DIT(6)) 

RETURN 

END 
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FUNCTION TEMP (ALT, TX, Tl, T3, T4, A2) 

C Function subprogram TEMP' calculates the temperature at altitude ALT 
C using equation (10) for altitudes between 90 and 125 km and equation 
C (13) for altitudes greater than 125 km, from SAO Report 313. 

C 

C Written by Mike Hickey, USRA 

PARAMETER BB = 4.5E-6 

U = ALT - 125. 

IF (U GT . 0.) THEN 

TEMP = TX + A2*ATAN(T1*U*(1. + BB*(U**2.5)) / A2) 
ELSE 

TEMP = TX + T1*U + T3*(U3**3) + T4*(U4**4) 

END IF 

END 


REAL FUNCTION MOL_WT*4 (A) 

C Subroutine MOL_WT' calculates the molecular weight for altitudes 
C between 90 and 105 km according to equation (1) of SAO report 313. 

C Otherwise, MOL_WT is set to unity. 

C 

C Written by Mike Hickey, USRA 
DIMENSION B (7) 

DATA B/ 28. 15204, -0.085586, 1.284E-4, -1.0056E-5, -1.021E-5, 
1.5044E-6, 9.9826E-8 / 

IF (A. GT. 105.) THEN 
MOLWT = 1 . 

ELSE 

U = A - 100. 

MOL_WT = B (1) 

DO 1 I = 2, 7 

MOL_WT = MOL_WT + B (I)*U**(I-1) 

1 CONTINUE 
END IF 

END 
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SUBROUTINE GAUSS (Zl, Z2, NMIN, R, TX, Tl, T3, T4, A2) 

C Subdivide total integration-altitude range into intervals suitable for 
C applying Gaussian Quadrature, set the number of points for integration 
C for each sub-interval, and then perform Gaussian Quadrature. 

C Written by Mike Hickey, USRA, NASA/MSFC, ED44, July 1988. 

REALM ALTMIN (9), C(8,6), X(8,6), MOL_WT 
INTEGER NG (8), NGAUSS, NMIN, J 

GRAVITY (ALTITUDE) = 9.80665 / (( 1. + ALTITUDE / 6.356766E3)**2) 

DATA ALTMIN/ 90., 105., 125., 160., 200., 300., 500., 1500., 2500. / 
DATA NG / 4, 5,6,6,6,6,6,67 


C Coefficients for Gaussian Quadrature ... 

DATA C / .5555556 , .8888889 , .5555556 , .0000000 , ! n=3 
.0000000 , .0000000 , .0000000 , .0000000 , ! n=3 
.3478548 , .6521452 , .6521452 , .3478548 , ! n=4 
.0000000 , .0000000 , .0000000 , .0000000 , ! n=4 
.2369269 , .4786287 , .5688889 , .4786287 , ! n=5 
.2369269 , .0000000 , .0000000 , .0000000 , ! n=5 
.1713245 , .3607616 , .4679139 , .4679139 , ! n=6 
.3607616 , .1713245 , .0000000 , .0000000 , ! n=6 
.1294850 , .2797054 , .3818301 , .4179592 , ! n=7 
.3818301 , .2797054 , .1294850 , .0000000 , ! n=7 
.1012285 , .2223810 , .3137067 , .3626838 , ! n=8 
.3626838 , .3137067 , .2223810 , .1012285 / ! n=8 


C Abscissas for Gaussian Quadrature ... 

DATA X/ -.7745967, .0000000, .7745967, .0000000 ,! n=3 
.0000000 , .0000000 , .0000000 , .0000000 , ! n=3 
-.8611363 ,-.3399810, .3399810, .861 1363 ,! n=4 
.0000000 , .0000000 , .0000000 , .0000000 , ! n=4 
-.9061798 , -.5384693 , .0000000 , .5384693 , ! n=5 
.9061798 , .0000000 , .0000000 , .0000000 , ! n=5 
-.9324695 , -.6612094 , - 2386192 , .2386192 , ! n=6 
.6612094 , .9324695 , .0000000 , .0000000 , ! n=6 
-.9491079 , -.7415312 , -.4058452 , .0000000 , ! n=7 
.4058452, .7415312, .9491079, .0000000 ,! n=7 
-.9602899 , -.7966665 , - 5255324 , -.1834346 , ! n=8 
.1834346 , .5255324 , .7966665 , .9602899 / ! n=8 
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R - 0.0 


DO 2 K = NMIN, 8 
NGAUSS = NG (K) 

A = ALTMIN (K) 

D = AMIN1 ( 72 , ALTMIN (K+l)) 

RR =0.0 
DEL = 0.5*(D - A) 

J = NGAUSS -2 
DO 1 1=1, NGAUSS 
Z = DEL*(X(I,J) + 1.) + A 

RR = RR + C(I,J)*MOL_WT(Z)*GRAVITY(Z) / TEMP (Z,TX,T1,T3,T4,A2) 

1 CONTINUE 
RR = DEL*RR 
R =R + RR 

IF (D EQ. Z2) RETURN 

2 CONTINUE 

RETURN 

END 
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c...test routine for j70met 
IMPLICIT NONE 

REAL TE, Z, TZ, DENS, DL 

TYPE *, 'Input exospheric temperature and altitude (km)' 
ACCEPT *, TE, Z 

CALL J70MET (TE, Z, TZ, DENS, DL) 

TYPE ' T = ' , TZ 

TYPE *, ' Density = ' , DENS 

TYPE log (density) = ' , DL 

STOP 

END 
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